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1.  Summary  of  ADCIRC,  v45.11,  Features 


The  Advanced  CIRCulation  Model  has  seen  tremendous  evolution  since  the  last 
transitioned  version,  v36.01.  The  code  has  been  translated  into  the  FORTRAN  90  language 
standard  and  the  two-dimensional  and  three-dimensional  codes,  once  separate,  are  now 
contained  within  a  single  unified  code.  The  2D  and  3D  unified  code  is  completely 
parallelized  using  the  MPI  protocol  and  metis  domain  decomposition.  The  code  structure 
itself  is  modular  with  individual  subroutines  containing  global  dimensioning,  cold  and 
hotstart  initialization,  parameter  specification,  wind  forcing,  harmonic  analysis,  three- 
dimensional  calculations,  and  the  main  time  stepping  loop.  Furthermore,  the  code  has 
become  multi-algorithmic  with  options  for  selecting  various  implementations  of  the 
conservative  or  non-conservative  Generalized  Wave  Continuity  Equation  (GWCE)  and/or 
momentum  equation  formulations. 

Other  notable  additions  to  the  ADCIRC  v45.11  model  code  are  options  to  use  a 
predictor-corrector  time  stepping  algorithm,  a  spatially  varying  GWCE  parameter  (TauO),  a 
nonlinear  hybrid  bottom  friction  relationship  and  a  wetting  and  drying  algorithm  has  been 
significantly  improved.  Under  appropriate  dynamical  conditions  the  predictor-corrector  time 
stepping  option  allows  time  steps  up  to  ten  times  larger  than  the  default  explicit  time 
stepping,  thereby  decreasing  the  computational  cost  of  a  particular  simulation.  Significant 
improvements  to  the  wetting/drying  capability,  including  the  option  of  a  nonlinear  hybrid 
bottom  friction  relationship,  have  been  implemented  to  improve  realism,  performance  and 
stability  of  the  model  when  wetting  and  drying  are  important.  Options  to  specify  a  spatially 
varying  GWCE  parameter  (TauO)  have  resulted  in  better  mass  conservation  and  stability  for 
highly  advective  flows. 

A  complete  description  of  the  theoretical  basis  for  ADCIRC  v44.xx  and  above  is 
given  in  http://adcirc.org/adcirc  theory  2004  12  08.pdf.  The  online  User’s  manual  for 
v45.11  is  found  at  http://adcirc.org/document/ADCIRC  title  page.html.  Below  are  detailed 
descriptions  of  the  latest  ADCIRC  code  features  including  wetting  and  drying,  a  hybrid 
nonlinear  bottom  friction  relationship,  predictor-corrector  time  stepping  and  the  spatially 
variable  GWCE  parameter  specification. 

a.  The  Wetting  and  Drying  Algorithm 

The  wetting  and  drying  scheme  is  based  on  a  simplified  one-dimensional  momentum 
balance  between  gravity  and  pressure  which  accounts  for  the  physics  of  the  flow,  along  with 
some  empirical  rules  to  help  ensure  stability.  Based  on  some  predetermined  criteria  (to  be 
discussed  later),  the  wetting  and  drying  ultimately  involves  classifying  elements  as  either 
wet,  which  means  there  is  a  sufficient  amount  of  water  over  that  element,  or  dry  meaning  that 
element  does  not  have  enough  water  over  it  to  be  a  part  of  the  active  computational  area. 

When  wetting  and  drying  of  elements  is  to  be  used,  ADCIRC  initializes  the  elevation 
at  all  nodes  such  that  the  total  water  depth,  Hj  at  node  j,  is  greater  than  or  equal  to  H0,  the 
minimum  specified  water  depth.  The  total  water  depth,  Hj  at  node  j,  is  defined  to  be  the 
bathymetric  depth,  z7  relative  to  mean  sea  level  (MSL)  plus  the  sea  surface  elevation,  /;,• 
relative  to  MSL,  see  f  igure  1.1. 
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Figure  1.1  Profile  of  bathymetry,  z,  and  free  surface  elevation,  r\,  relative  to  mean  sea  level 
(MSL). 


Figure  1.2.  Wet/Dry  interface  acts  as  a  no  normal  flow  barrier.  Red  dots  marked  with  a  "D” 
indicate  dry  nodes  and  red  dots  marked  with  a  "W”  indicate  wet  nodes. 

The  wetting  and  drying  algorithm  can  be  broken  down  into  two  main  components,  the 
first  is  a  drying  phase  and  the  second  is  a  wetting  phase.  The  drying  phase  is  based  on  nodal 
criteria.  During  this  phase,  the  total  water  height,  H,  at  all  active  nodes  are  examined  to 
determine  if  there  is  enough  water  present  for  them  to  stay  active,  if  not,  those  nodes  are 
dried,  see  figure  1.3a  The  threshold  value  between  a  classification  of  wet  or  dry,  is  based  on 
a  total  water  depth  greater  than  or  equal  to  Hq.  If  the  total  water  depth  is  at  or  below  Ho  then 
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the  node  is  classified  as  dry.  Furthermore,  if  the  total  water  depth  at  the  now  dry  node  has 
fallen  below  an  absolute  minimum  water  depth,  Habsm.  =  0.8 Ho,  then  elevations  are  reset  such 
that  the  total  water  height  is  Habsm-  Thus  newly  dried  nodes  will  have  a  total  water  depth  that 
is  bounded  below  by  Habsm  and  above  by  Ho.  This  variance  in  water  height  can  explain  some 
of  the  asymmetry  observed  during  re-wetting  phases  along  a  front.  This  will  be  elaborated 
on  in  the  discussion  of  results  (Section  3). 

The  second  phase  of  the  wetting  and  drying  algorithm  is  the  wetting  phase,  see  figure 
1.3b  for  a  flow  chart  illustration.  It  is  in  this  phase  that  we  evaluate  the  simplified 
momentum  balance  between  gravity  and  pressure, 


where  g  is  acceleration  due  to  gravity.  Upon  rearrangement  of  (1)  we  get 

~  (  \ 

g 


v  = 


drj_ 

dx 


Cn 


where  r,  =  — .  We  will  refer  to  r,  as  a  bottom  stress  coefficient. 

*  H  k 


(2) 


In  the  wetting  phase  of  the  wetting  and  drying  algorithm,  instead  of  individual  nodes 
being  considered,  dry  elements  are  examined.  An  element  is  considered  wet  only  if  all  its 
nodes  are  wet.  All  dry  elements  that  have  exactly  two  wet  nodes  are  examined  to  determine 
if  they  should  become  wet  and  thus  force  their  one  dry  node  to  become  wet.  Let’s  assume 
that  we  have  an  element  with  two  wet  nodes,  label  them  1  and  2  and  only  one  dry  node,  label 
it  k.  Before  wetting  can  occur,  the  total  water  column  height  at  the  two  wet  nodes  must  be 
greater  than  or  equal  to  Hofj.  =  1 2H0  .Using  the  higher  value  for  wetting  than  for  drying 

reduces  the  chances  of  re-wetting  nodes  that  have  just  dried,  which  could  result  in 
instabilities  due  to  nodes  turning  on  and  off.  Let  us  assume  that  both  of  the  wet  nodes  meet 
this  criterion.  Now  the  simplified  momentum  balance  is  employed  to  determine  if  there  is 
enough  water  to  wet  the  element.  First  the  wet  node  with  the  largest  elevation  will  be  used  in 
the  computation,  without  loss  of  generality  assume  it  is  local  node  number  one,  Hi,  The 
velocity  that  would  result  in  the  water  moving  from  node  1  to  node  k  is  computed  using  a 
discrete  form  of  equation  (2), 


^wet  g 


A  77 

vAxy 


kwet 


(3) 


with  A  77  =(r/i~rik).  Ax  the  distance  between  node  1  and  k  and  a  bottom  stress  coefficient, 
Tkwet  ’  computed  using  either  a  quadratic  or  hybrid  nonlinear  bottom  friction  relationship  for 


Cd  with  the  velocity  set  at  the  prescribed  minimal  velocity  allowed,  v  =  vmin ,  and  H  =  Hl. 


If  vwet  >  vmin  then  the  element  is  reclassified  as  wet  and  is  added  to  the  computational 

domain.  Notice  that  the  simplified  momentum  balance  does  not  take  into  account  any 
directionality  of  flow.  Thus,  one  could  (and  does)  have  cases  where  the  flow  is  being  forced 
toward  one  direction  and  wetting  is  occurring  in  another  direction(s).  We  will  see  evidence 
of  this  in  the  first  case  study  presented  in  Section  3. 
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There  are  a  few  additional  empirical  rules  that  make  up  the  wetting  and  drying 
algorithm  which  deal  mostly  with  weir  overtopping.  These  rules  are  not  effected  by  Ha  and 
v  ■  and  thus  will  not  be  discussed  here. 


Figure  1 .3.  Decision  tree  for  (a)  drying  a  wet  node  and  (b)  wetting  a  dry  element. 

When  dry  elements  are  made  wet,  the  water  required  to  fill  that  element  is  just  added 
to  the  overall  water  mass  in  the  active  computational  area.  Likewise,  when  a  wet  element  is 
made  dry,  whatever  amount  of  water  was  left  in  the  now  dry  element  is  just  removed  from 
the  total  water  mass.  This  has  immediate  negative  consequences  when  coupling  ADCIRC 
solutions  with  wetting  and  drying  to  a  set  of  transport  equations.  Some  might  muse  that 
since  many  of  the  domains  have  an  open  boundary,  e.g.,  the  open  ocean  or  a  river,  that  is 
actively  forcing  the  simulation,  then  how  important  is  accounting  for  mass  balance?  The 
answer  depends  on  the  application.  If  one  primarily  interested  in  high  water  mark  values 
from  a  hurricane  storm  surge,  such  mass  imbalances  would  not  be  important.  But  for 
modeling  small-scale  non-extreme  events,  such  as  tidal  flooding  and  drying,  the  lack  of  mass 
conservation  may  be  detrimental  especially  when  transport  processes  are  of  great  importance. 
Furthermore,  as  we  will  see  in  the  cases  studies  presented  in  Section  4  the  sudden  addition  or 
subtraction  of  water  mass  during  the  solution  process  can  produce  near  shocks  in  the  velocity 
solutions.  These  shocks  can  be  strong  enough  to  cause  model  instabilities  so  severe  as  to 
result  in  premature  model  termination. 

Practically,  wetting  and  drying  within  the  ADCIRC  model  is  invoked  with  a 
parameter  specification  of  NOLIFA,  the  parameter  controlling  the  nonlinear  finite  amplitude 
terms,  equal  to  2  or  3.  Under  wetting  and  drying  values  of  the  minimum  water  depth,  H0 , 
represent  the  nominal  water  depth  at  a  node  for  it  to  be  considered  dry,  typically  0.01  -  0.1 
m.  Additional  parameters,  NODEDRYMIN,  NODEWETMIN,  and  VELMIN,  are  also 
required  to  activate  wetting  and  drying.  See  the  online  user’s  manual  for  details  on  their 
specification  (http://adcirc.org/document/ADCIRC  title  page.html)  and  see  Section  3a  for  a 
discussion  of  the  sensitivity  of  these  parameters  on  the  performance  of  the  wetting  and  drying 
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algorithm. 

b.  Nonlinear  Bottom  Friction  Relationships 


Important  to  the  physical  wetting  and  drying  process  is  the  discrete  representation  of 
bottom  friction.  ADCIRC  has  an  option  for  three  types  of  bottom  friction, 

h=CDv  (4) 

defined  by  the  choice  of  the  drag  coefficient  Cd .  The  first  choice  for  bottom  friction 
coefficients  is  a  linear  friction  relationship  CD  =  cf ,  the  second  one  is  a  quadratic  friction 

relationship, 

cD=cfU  (5) 


and  the  third  is  a  hybrid  nonlinear  bottom  friction  which  results  from  using  a  depth 
dependent  drag  coefficient, 


CD  =  max 


C 


D  min 


l+ 


H, 


break 


H 


,10 


,-4 


(6) 


where  Comm  is  either  of  the  previously  defined  linear  or  quadratic  friction  drag  coefficients,  v 
is  the  water  velocity,  and  H  is  the  total  water  column  depth.  The  depth  dependent  drag 
coefficient  was  introduced  into  ADCIRC  in  the  report  by  Luettich  and  Westerink  (1995)  and 
is  discussed  in  detail  in  that  report.  In  brief,  when  in  deep  water,  H  >  Hbreak ,  the  drag 

coefficient  approaches  Comm,  the  standard  linear  or  quadratic  drag  coefficient  relationships. 
However,  when  in  shallow  water,  H  <  Hbreak ,  Cd  approaches  CDmm  (Hbreak  / H\  .  According  to 

Luettich  and  Westerink  (1995)  the  exponent  0  determines  how  rapidly  CD  approaches  each 
asymptotic  limit,  while  the  exponent  y  determines  the  rate  at  which  the  friction  coefficient 
increases  as  the  water  depth  decreases.  Using  the  hybrid  nonlinear  bottom  friction  case  when 
wetting  and  drying  means  that  the  friction  coefficients  will  be  higher  than  those  computed 
when  using  the  standard  linear  or  quadratic  option.  The  suggested  default  values  for 
ADCIRC's  hybrid  nonlinear  bottom  friction  are  Hbreak=2m,  0  =10,  y  =  1.3333,  and  Cf  = 
0.0030. 

The  choice  of  nonlinear  bottom  friction  is  determined  by  the  value  of  NOLIBF  in  the 
fort.  15  parameter  file.  NOLIBF  =  1  results  in  the  standard  nonlinear  quadratic  friction  law, 
equation  (5),  while  NOLIBF  =  2  invokes  the  hybrid  nonlinear  bottom  friction  formulation  of 
equation  (6). 


c.  Predictor-Corrector  Time  Stepping 


Currently,  non-linear  applications  with  ADCIRC  have  stability  issues  unless  a  severe 
Courant  number  restriction  is  imposed.  The  Courant  number  is  defined  as 

c-xr  <7) 


where  c  =  ,Jgh  is  the  linear  wave  celerity,  Ax  is  minimum  node  spacing  and  At  is  the  time 

step.  In  practice,  we  have  found  that  for  deep  ocean  flows,  a  practical  upper  bound  of  the 
Courant  number  (C,  )  is  0.5  in  order  to  maintain  stability;  however,  an  even  tighter  constraint 
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(e.g.  Cr  «  0.1 )  must  be  imposed  if  the  simulation  includes  barrier  islands,  constricted  inlets, 
or  wetting  and  drying  of  near-shore  elements.  In  order  to  relax  this  restriction,  an  alternative 
time-marching  procedure  was  proposed  that  treats  the  non-linear  terms  implicitly  (Kolar  et 
al.,  1998). 

Within  the  standard  ADCIRC  model  semi-implicit  time  marching  algorithm,  linear 
terms  are  evaluated  implicitly  and  the  non-linear  terms  explicitly.  At  the  past  and  present 
time  levels  in  ADCIRC,  elevation  and  velocity  values  are  known  (either  from  initial 
conditions  or  previous  calculations).  The  original  algorithm  takes  the  elevation  and  velocity 
values  for  the  past  (k  -1)  and  the  present  ( k )  and  uses  them  to  calculate  the  values  for  the 
future  ( k+1 )  time  level  for  the  linear  terms.  However,  the  non-linear  terms  are  evaluated 
using  only  the  elevation  and  velocity  values  at  the  present  time  level  ( k ).  Kolar  et  al.  (1998) 
hypothesized  that  the  stability  constraint  stems  primarily  from  this  explicit  evaluation  of  non¬ 
linear  terms.  A  predictor-corrector  time-marching  algorithm  is  introduced  to  evaluate  the 
non-linear  terms  implicitly.  The  predictor  stage,  which  is  equivalent  to  the  original  algorithm, 
evaluates  the  non-linear  terms  using  values  from  the  present.  Predicted  future  values,  called 
k *,  and  the  already-known  present  (k)  and  past  (k~  1)  values  are  then  used  to  obtain  corrected 
values  for  the  future  (k  +  1)  time  level.  The  corrector  stage  can  be  repeated  as  many  times  as 
necessary  until  convergence.  In  all  applications  to  date,  a  single  iteration  of  the  corrector 
stage  appears  to  be  sufficient. 

Through  the  use  of  time  weight  coefficients,  users  have  the  option  to  distribute  the 
relative  contribution  of  the  non-linear  terms  over  the  three  time  levels.  Comprehensive  ID 
studies  by  Dresback  and  Kolar  (2001)  and  accompanying  2D  studies  (Dresback  et  al.  2004) 
have  shown  that  optimal  coefficients  are  problem  dependent,  but  that  near-optimal  results  for 
any  domain  are  found  by  centering  the  GWCE  time  weights  at  k  (meaning  that  the  time 
weights  for  the  non-linear  terms  are  weighted  equally  between  k  +1  (or  k*),  k,  k  - 1 )  and 
centering  the  non-  conservative  momentum  time  weights  at  k  +  14  (meaning  that  the  terms 
are  weighted  equally  between  k  and  k  +1  (or  k *)).  Practically,  the  predictor-corrector  time 
marching  is  activated  through  the  use  of  a  negative  value  for  the  time  step,  DTDP. 

d.  Spatially  Variable  GWCE  Parameter,  TauO 

The  Generalized  Wave-Continuity  equation  (GWCE)  formulation  used  in  the 
ADCIRC  model  contains  the  parameter,  t0  or  TauO,  that  serves  as  a  weighting  between  the 
primitive  continuity  and  the  wave  portions  of  the  GWCE  equation  (Kinnmark,  1996).  A 
value  for  r0  equal  to  zero  results  in  a  pure  wave  equation  while  a  value  of  1  for  ra  leads  to 
the  primitive  continuity  equation.  The  value  selected  for  t0  strikes  a  balance  between 
stability  and  mass  conservation.  A  pure  primitive  continuity  formation  (  t0  =  1 )  results  in  an 
unstable  solution  whereas  tq  =  0  leads  to  the  greatest  mass  imbalance.  ADCIRC  now 
provides  an  option  to  specify  spatially  variable  values  of  this  weighting  coefficient  so  as  to 
minimize  mass  conservation  issues  in  shallow  water  regions  or  in  regions  where  wetting  and 
drying  is  likely  to  occur. 

Spatially  variable  values  of  t0  are  automatically  determined  and  invoked  using  a 
flagged  value  for  t0  that  is  negative.  Under  the  spatially  variable  za  option,  the  value  is 
specified  as  follows: 

Depth  >  10.0,  r0  =0.005 
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Depth  <  10.0,  t0  =0.020 


This  approach  leads  to  smaller  mass  errors  in  the  shallower  waters  while  maintaining 
stability  over  deeper  waters. 

2.  Testing  Strategy 

Objectives  of  an  ADCIRC,  v45.11,  model  application  are  to  represent  two- 
dimensional  and  three-dimensional  currents,  water  levels  and  inundation  extent  under 
barotropic  dynamical  conditions  in  coastal  regions.  Barotropic  dynamics  are  present  when 
the  density-driven  component  of  the  flow  is  minimal  or  non-dominant.  There  are  many 
circumstances  in  the  coastal  ocean  wherein  barotropic  conditions  prevail. 

The  validation  test  cases  selected  aim  to  span  the  range  of  possible  applications  of  the 
ADCIRC  model  under  these  circumstances  in  order  to  evaluate  the  model’s  ability  to  predict 
realistic  and  accurate  coastal  dynamics.  Furthermore,  many  of  the  test  cases  are  identified 
with  established  or  published  benchmark  test  cases.  In  most  instances  field  measurements  are 
available  and  descriptions  of  the  test  cases  can  be  found  in  the  literature.  Table  1  presents  the 
matrix  of  validation  test  cases  and  the  component  of  the  barotropic  flow  under  evaluation. 


Table  1 .  Description  of  the  Validation  Test  Cases  for  ADCIRC  v45. 1 1 . 


Test  Case 

Testing  Parameters 

Physics 

Elevation 

Wetting/Drying 

Currents 

Dimensionality 

Idealized 

Beaches 

inundation 

V 

V 

2D 

Hurricane 

Katrina 

storm  surge, 
inundation 

V 

V 

V 

2D 

Delaware  Bay 

tides  and 
winds 

V 

V 

2D  and  3D 

Rattray  Island 

tides 

V 

3D 

MREA07 

Software 

V 

V 

2D  and  3D 
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3.  Validation  Test  Results  -  Wetting  and  Drying 

a.  Parameter  Sensitivity  Analysis 

To  gain  a  greater  understanding  of  how  the  control  parameters  H0  and  Vmin  effect 

the  wetting  and  drying  process,  the  simplified  momentum,  equation  (3),  is  analyzed  through  a 
series  of  experiments. 

Recall  that  the  discrete  wetting  velocity,  given  in  equation  (3),  must  be  greater  than  or 
equal  to  the  minimum  velocity  set  by  the  user,  vwet  >  Vmin ,  in  order  for  wetting  to  occur.  For 

simplification,  the  standard  nonlinear  bottom  friction  relation  (4)  is  assumed.  This 
assumption  will  not  alter  the  main  analysis  shown  below.  To  determine  the  relational  effects 
that  Vmin  has  on  Vwet  we  differentiate  equation  (3)  with  respect  to  Vmin , 


dVwet  J 

'Ari' 

-1  1 

_  _  ^wet 

dv  .  ~8\ 

min 

vAry 

v  •  r, 

y  min  kwet  ) 

V 

mm 

Since,  vwet  >  0  and  Vmin  >  0 ,  the  negative  sign  on  the  derivative  implies  an  inverse 

relationship  between  Vmin  and  Vwet ,  namely,  decreasing  the  value  of  Vmin  would  have  the 

effect  of  increasing  the  computed  wetting  velocity.  The  result  is  two  paths  to  easier  wetting 
of  the  element.  First  we  have  a  smaller  velocity  threshold  to  overcome  with  a  decreased  value 
of  Vmin .  Secondly,  the  decrease  in  Vmin  also  results  in  a  higher  computed  velocity,  as  a 

result  of  the  decreased  drag  coefficient  from  f . 

In  order  to  determine  the  effects  H0  has  on  Vwet ,  we  will  differentiate  Vwet  with 
respect  to  H0 .  But  first  Vwet  must  be  expressed  in  terms  of  H0 .  Assume  that  we  have  a 
dry  element  that  has  two  wet  nodes  and  one  dry  node  that  has  never  been  wet.  This  is  to 
ensure  that  its  total  water  depth  is  H0 .  Without  loss  of  generality,  assume  that  local  node 
number  1  is  a  wet  node  with  the  largest  elevation  and  the  dry  node  is  local  node  number  k. 
Since  node  1  is  wet  we  know  that  Hl  >  H ahsm  =  0.8//0  and  in  order  for  node  k  to  be 

considered  for  wetting,  H]  >  H off  =  1 .2 H0 .  Furthermore,  since  node  k  is  dry  and  has  never 
been  wet,  it's  elevation  has  been  set  such  that  Hk  =  H0  therefore  rjk=H0-zk.  In  order  for 
wetting  to  occur,  A.rj  =  rjl—T]k  must  be  greater  than  0,  otherwise  vwet  =0.  The  elevation  at 
local  node  1  can  be  written  as  zk  and  Hl  can  be  written  as  Hl=H0+  SM ,  where 

8h j  is  a  function  of  elevation,  l) .  In  fact,  since  H]  >  1 .2H(j  we  know  that  <5^  >  0.2Hf) . 
Combining  these  we  can  write  r/x  =  H0+Shl  —zx ,  which  we  use  in  decomposing  Arj ,  namely, 

A?7  =  Vi  ~Vk  =  8h\  +(zk  ~zi)  (9) 

Again  for  the  sake  of  convenience,  a  nonlinear  bottom  friction  term  is  assumed  instead  of  the 
hybrid  nonlinear  one.  The  discrete  wetting  velocity,  equation  (3)  as 
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^ wet  § 
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1  kwet 


=_g_ 
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(Vx-Vu) 
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CrV  ■ 

V  /  min 
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^C/Vm  in 
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(10) 


Now  in  order  to  determine  the  relational  effects  of  H0  and  Vwet ,  we  differentiate  equation 
(10)  with  respect  to  H0  to  obtain 


<9v, 


wet  — 


g 


dHn 


Axc,v  . 

J  min 


{sh\+(zk~z\)) 


(ii) 


Now,  by  assumption,  A//  >  0  so  that  Shl  +(zk  —  z,)  >0 .  Therefore  the  nonnegative 
derivative  in  equation  (11)  indicates  a  direct  relationship  between  H0  and  Vwet ,  namely 
increases  in  H0  mean  increases  in  Vwet .  Using  the  hybrid  nonlinear  bottom  friction  only 
alters  the  magnitude  of  the  relationship,  not  the  sign  of  the  derivative.  Overall  an  increase  in 
H0  means  that  initially  it  is  harder  to  wet  an  element,  however,  once  the  wetting  begins, 
further  wetting  should  be  easier  to  obtain.  A  larger  H0  requires  water  to  pile  up  higher 
along  the  wet  edge  of  the  dry  element,  see  figure  1.2.  The  resulting  higher  water  level 
increases  the  forward  velocity  of  the  water  because  the  change  in  elevation,  A  rj  expressed  in 
equation  (9)  would  also  increase  because  Sh]  >  0.2Htl . 


The  effects  of  H0  and  Vmin  on  the  WAD  algorithm's  performance,  can  also  be 

examined  by  considering  the  behavior  of  Vwet  in  terms  of  A  77.  We  will  graphically  show 

the  effects  that  H0  and  Vminhave  in  terms  of  a  required  change  in  elevation,  A  77,  that  is 

large  enough  to  enable  wetting  to  occur,  i.e.  so  that  vwet  >  vmin .  This  analysis  is  important 

when  considering  how  to  select  appropriate  values  for  the  parameters  H0  and  Vmin  given  a 

particular  mesh  and  topography.  The  bathymetric  gradients  shown  in  equation  (9)  can  be 
figured  out  in  advance  and  the  user  can  get  an  estimate  of  how  much  water  has  to  build  up 
before  wetting  can  occur  for  given  parameter  choices. 

For  this  analysis  we  assume  that  the  total  water  depth  at  the  wet  neighbor  node, 
Hl  =1.2 H0,  is  the  minimal  allowable  value  for  wetting  to  occur.  For  a  mesh  with  a  nodal 

spacing  of  500m  we  show  in  figure  2.1a  the  contours  of  the  required  elevation  changes  for 
both  the  nonlinear  and  the  hybrid  nonlinear  bottom  friction  cases.  Notice  that  the  slopes  of 
the  contours  in  figure  2.1a  are  both  positive  with  respect  tovmin ,  however  their  curvature  is 

different.  The  nonlinear  bottom  friction  case  has  a  positive  curvature  while  the  hybrid 
nonlinear  bottom  friction  case  has  a  negative  curvature.  This  is  pointed  out  to  support  our 
rational  for  performing  the  analysis  using  the  less  complicated  nonlinear  bottom  friction  case. 
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Furthermore,  the  nonlinear  bottom  friction  case  allows  for  smaller  Rvalues,  hence 

smaller  H0,  than  the  hybrid  nonlinear  bottom  friction  case  for  the  same  required  change  in 
elevation.  Thus  as  expected,  the  nonlinear  bottom  friction  case  allows  wetting  to  occur  more 
easily.  Mesh  element  size  acts  as  a  linear  multiplier  on  the  elevation  change  required  for 
wetting  to  occur.  The  smaller  the  mesh  size  the  easier  wetting  can  occur.  This  can  be  clearly 
observed  in  figure  2.1b  where  for  the  hybrid  nonlinear  bottom  friction  case,  the  effects  that 
H0  and  vminhave  on  a  required  elevation  change  of  0.1m  are  shown  for  meshes  with 

element  sizes  of  Ax  =  500,  250,  and  125  m. 
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Figure  2.1:  Wetting  and  drying  parameter  sensitivity  analyses  comparing  (a)  the  required 
change  in  elevation  relative  to  Vmin  and  H0  using  nonlinear  bottom  friction  and  hybrid 

nonlinear  bottom  friction  relationships  and  (b)  the  effect  that  mesh  spacing  has  on  the 
required  elevation  change. 

As  was  shown  by  Dietrich  et  al.  (2004,  2006),  an  increasing  value  of  the  free 
parameter  H0  also  increases  the  mass  balance  errors.  The  impact  that  H0  has  on  global 
mass  balance  is  addressed  via  a  test  case  in  two  dimensions,  the  Enclosed  Circular  Basin , 
presented  in  Section  3.2. 


For  all  wetting  and  drying  test  cases,  the  2D  depth  integrated  formulation  of 
ADCIRC,  version  45.11,  is  applied.  The  pertinent  input  parameter  values  (fort.  15)  are 
specified  below.  If  the  input  values  deviate  from  those  listed  here,  the  change  will  be  noted 
under  the  description  of  each  test  problem.  The  ADCIRC  model  input  file  (fort.  15)  parameter 
specifications  are: 
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IM  =  0  -  Model  selection  parameter 

NOLIBF  =  2  -  Bottom  friction  term  selection  parameter 

NOLIFA  =  2  -  Finite  amplitude  term  selection  parameter 

NOLICA  =  0  -  Spatial  derivative  convective  selection  parameter 

NOLICAT  =  0  -  Time  derivative  convective  term  selection  parameter 

NWP  =  0  -  Variable  bottom  friction  and  lateral  viscosity  option  parameter 

NCOR  =  0  -  Variable  Coriolis  in  space  option  parameter 

NTIP  =  0  -  Tidal  potential  option  parameter 

NWS  =  2  -  Wind  stress  and  barometric  pressure  option  parameter 

NRAMP  =  1  -  Hyperbolic  Tangent  Ramp  Function  Option 

TAUO  =  -0.020  -  Weighting  factor  in  GWCE 

0.35  0.30  0.35  =  Time  weighting  factors  for  the  GWCE  equation 

FF ACTOR, HBREAK,FTHETA,FGAMMA  =  0.0030  2.0  10  1.33333 

ESL  =  0.50  -  Lateral  eddy  viscosity  coefficient 

CORI  =  0.0  -  Coriolis  parameter 

Setting  IM  =  0  enables  a  barotropic  2D  depth  integrated  run  using  the  original  GWCE  and 
momentum  equation  formulations.  Setting  NOLIBF  =  2  enables  the  hybrid  nonlinear  bottom 
friction  calculation.  Setting  NOLIFA  =  2  enables  the  wetting  and  drying  algorithm. 
Specifying  NRAMP  =  1  enables  a  hyperbolic  tangent  ramping  function,  see  equation  (12) 
that  is  used  to  ramp  up  the  input  values  from  an  initial  factor  of  zero  to  nearly  96  percent  of 
the  specified  input  by  the  end  of  the  ramp  duration,  given  by  the  input  variable  DRAMP. 

RAMP  =  tanh^.O*  IT  *  DT  /  (865400.*  DRAMP))  j  (12) 


where  IT*DT  is  the  elapsed  time  in  seconds  of  the  simulation,  and  DRAMP  is  the  effective 
duration  of  the  ramp  given  in  days. 

A  complete  description  of  these  and  all  other  ADCIRC  input  variables  is  given  on  the 
ADCIRC  website,  http://www.adcirc.org. 

b.  Enclosed  Circular  Basin 

The  enclosed  circular  basin  of  Xie  et  al.  (2004)  is  used  to  examine  the  non¬ 
conservative  nature  of  the  ADCIRC  wetting  and  drying  algorithm.  The  physical  domain  is  a 
120A:m  by  120 km  square  with  a  circular  basin  of  radius  30 km  located  at  the  center,  see  figure 
2.2a.  The  domain  is  discretized  using  structured  right  triangular  elements  formed  by 
diagonally  cutting  uniform  squares  of  size  600  meters.  The  bathymetry/topography  set  up  is 
the  same  as  in  Xie  et  al.  (2004)  where  the  depth  at  the  center  of  the  bowl  is  9.0  meters  and 
decreases  linearly  to  1.0  meter  at  the  land/sea  interface  along  the  rim  of  the  bowl.  The  land 
elevations  begin  with  a  value  of  -0.25  meters  at  the  rim  and  change  linearly  with  a  constant 
slope  of  -1/7500  to  the  exterior  boundaries.  The  outer  most  boundaries  were  set  to  a 
mainland  boundary  type  with  no  normal  flow  as  an  essential  boundary  condition  and  free 
tangential  slip  allowed.  The  same  wind  forcing  as  Xie  et  al.  (2004)  derived  from  a  Holland 
(1980)  axis-symmetric  hurricane  model  were  used  to  force  the  run.  Note  a  correction  to 
equation  11  in  Xie  et  al.  (2004)  where  the  wind  drag  Cd  coefficients  should  have  been 
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10  3c  =0.62  +  1.56 


K 


for  wind  speeds  1  < 


K 


<  10  m/s. 


Wind  fields  were  supplied  every 


30  minutes  to  ADCIRC  for  a  period  of  54  hours,  see  figure  2.2b  for  storm  center  locations 
relative  to  the  computational  domain.  The  simulation  was  for  54  hours  using  a  6  hour  ramp 
period  for  the  wind  forcing  and  a  time  step  of  5.0  seconds.  A  series  of  experiments  is 
performed  in  which  the  wetting  and  drying  control  parameters  H0  and  V  •  are  varied.  The 


values  presented  are  for  H0  =  (0.1,  1.0,  5.0)  meters  and  for  Vmin  =0.01,  0.05,  0.1,  0.2,  0.5, 
1.0)  meters  per  second. 
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Figure  2.2.  The  (a)  topography/bathymetry  of  the  enclosed  circular  basin  and  (b)  the  storm 
track  locations,  in  red,  relative  to  the  computation  domain  shown  in  blue. 


In  figure  2.3,  a  series  of  contour  plots  showing  the  wet/dry  interface  for  various 
values  of  V  .  at  times  18  hours,  24  hours,  and  54  hours  into  the  simulation  (rows)  and 

various  H0  values  (columns).  In  each  of  these  images  the  original  coastline  of  the  basin  is 
drawn  as  a  frame  of  reference  from  which  to  view  the  progress  of  the  wetting  and  drying 
front.  The  first  two  rows  show  the  initial  wetting  phase  of  the  simulation  while  the  last 
image  shows  what  should  be  a  drying  phase  with  all  the  water  returning  to  the  original  basin 
under  minimal  winds.  Examining  the  images  at  the  18  hour  mark,  we  see  that  by  increasing 
H0  the  drying  front,  along  the  northeast  border  of  the  basin,  moves  further  toward  the  bowl's 
center  away  from  the  original  coastline.  This  is  to  be  expected,  due  to  the  drying  criteria  of 
the  WAD  algorithm.  The  effects  of  H0  on  the  wetting  front,  along  the  southwest  portion  of 
the  basin,  are  more  complicated  to  analyze.  For  the  smallest  H0 ,  the  wetting  front  is  more 
diffuse  and  less  advanced  away  from  the  bowl's  center  than  the  H0  =  0.5 m  case.  We 
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attribute  this  behavior  to  the  fact  that  the  smaller  H0  allows  the  water  to  meander  around  the 
concentric  rings  of  bathymetry  more  than  higher  values  of  H0  whose  wetting  fronts  are 
more  in  line  with  the  direction  of  the  wind  forcing.  By  comparing  the  contours  at  the  1 8  hour 
mark,  for  H0  =  0.1  m  and  H0  =  1.0 m,  we  notice  that  for  higher  H0,  the  wetting  front  takes 
longer  to  begin  advancing.  This  is  because  higher  H0  values  require  more  water  to  build  up 
behind  the  wetting  front  before  allowing  the  wetting  front  to  advance.  At  the  18  hour  mark 
for  the  cases  with  H0  =  1 .0 m ,  the  water  levels  are  still  building  and  have  not  started  to 

propagate  outside  the  basin.  By  examining  the  24  hour  mark  for  all  H0  values,  as  well  as 
the  18  hour  mark  for  H0  =0.1  m  and  H0  =0 .5m,  we  see  that  for  increasing  H0  values,  that 
once  the  water  begins  to  move,  the  wet/dry  front  moves  further  and  is  more  focused  in  the 
direction  towards  which  the  winds  are  actually  blowing.  At  the  54  hour  mark,  we  again  see 
that  for  increased  H0  values,  the  faster  the  drying  process  acts.  In  fact  for  the  smallest  Ha 
value,  there  is  still  a  large  area  of  land  that  is  inundated,  even  though  the  wind  forcing  is 
negligible.  The  case  for  H0  =1.0 m  and  Vmin  =0.2  m/s  at  the  54  hour  mark  seems  to  have 

settled  to  an  equilibrium  solution  outside  the  original  basin,  for  which  we  do  not  have  a  clear 
explanation. 

The  effects  of  Vmin  on  the  wetting  front  can  be  seen  clearly  in  the  plot  for  Ha  =  0. 1  m  at  24 
hours,  with  the  wetting  front  being  further  advanced  for  decreasing  values  of  Vmin  ,  as 
predicted  by  our  previous  analysis.  Again,  a  smaller  value  of  Vmin  results  in  a  larger,  more 

diffuse  wetting  front.  This  same  effect  is  observed  for  the  other  H0  values  at  24  hours  and 
for  the  1 8  hour  mark  with  Ha  =  0.  lm  and  H0  =  0. 5m . 

In  figure  2.4  are  shown  the  time  series  of  elevation  at  three  nodal  points,  whose 
locations  are  depicted  in  (a),  for  H0  =  0.5 m  and  1.0m  and  for  a  fixed  Vmjn  =  0.05  m/s.  When 

a  node  dries,  the  elevation  solution  is  represented  using  a  value  of  0.0m  to  better  distinguish 
between  wet  and  dry  locations.  Notice  that  for  all  locations  the  final  elevations  are  near 
H  ,  ,  which  means  the  final  water  elevations  are  higher  than  the  initial  water  elevation 

even  though  the  domain  is  closed  and  the  water  is  by  this  time  nearly  at  rest.  This  is  a  direct 
consequence  of  the  nodal  drying  criteria  in  the  WAD  algorithm.  Notice  that  for  node  1 
(figure  2.4b)  located  in  the  center  of  the  basin  in  the  deepest  water,  the  elevation  solution  is 
relatively  smooth,  with  oscillations  occurring  only  when  the  water  level  has  settled  in 
around//  ,  .  Also,  note  that  node  1  never  dries  and  that  the  elevation  solution  for  the  case 

with  Hn  =  1.0m  returns  to  H  ,  faster  than  in  the  case  when  Hn  =  0.5m .  At  node  numbers 

51  and  10252,  the  elevation  solutions  (figures  2.4c  and  2.4d)  both  experience  periods  of 
drying  and  then  re-wetting.  The  higher  H0  values  result  in  drying  occurring  faster,  and  in 
the  case  of  node  10252,  results  in  wet/dry  oscillation  during  the  first  6  hours  of  the 
simulation.  For  both  of  these  locations,  the  maximum  elevation  reached  during  a  wetting 
phase  is  not  affected  significantly  by  the  value  of  H0 ,  see  figure  2.4c  around  24  hours  and 
figure  2.4d  around  28  hours.  However,  during  the  final  drying  stage,  after  40  hours,  both  of 
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these  nodes  experience  larger  oscillations  in  elevation  around  H  ,  than  node  number  1 

during  the  same  time,  for  H0  =1.0 m.  These  oscillations  are  possibly  due  to  the  drying  of 
neighboring  nodes. 

Figure  2.5  contains  a  time  series  plot  of  the  total  number  of  wet  nodes  in  (a)  and  in  (b) 
the  corresponding  computed  total  water  volume  in  the  domain  for  both  values  of H0. 
Obviously,  the  WAD  is  not  mass  conserving  even  after  all  the  water  settles  back  into  the 
basin.  The  total  water  height  is  now  bounded  between  H absm  and  H0  resulting  in  a  higher 

resting  elevation.  Notice  that  for  the  H0  =  0.5 m  case  wetting  began  outside  the  basin  earlier, 
wetted  more  nodes,  and  took  longer  to  drain  back  into  the  basin  than  the  Ha  =1.0 m  case. 
This  is  again  due  to  the  wetting  criteria  in  the  WAD,  where  the  smaller  H0  allows  for  more 
meandering  along  the  concentric  bathymetry  rings. 
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Figure  2.3.  Contours  of  the  wetting  and  drying  front  at  18,  24  and  54  hours  (rows)  into  the 
simulation,  using  different  H0  values  (columns),  and  different  Vmin  values  represented  by 

different  color  contours  lines.  The  original  coastline  of  the  basin  (in  black)  is  included  as  a 
frame  of  reference. 
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Elevation  at  node  #1 


Time  (Hours) 


(c) 


Elevation  at  node  #51 


(d) 


Figure  2.4.  The  (a)  location  of  nodal  points  along  the  enclosed  circular  basin  with  their  time 
series  elevations  (b-d)  relative  to  the  geoid,  using  H0  =  0.5  and  H0=  1.0  with  V  .  =  0.05  . 


Figure  2.5.  In  (a)  the  number  of  wet  nodes  and  in  (b)  the  computed  total  volume  of  water 
during  the  simulation  for  two  Ha  values,  H0  =  0.5  and  H0  =  1.0 ,  with  V  .  =  0.05  . 
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c.  Idealized  Inlet 


The  impact  of  the  WAD  algorithm's  two  free  parameters  on  the  velocity  solution  for  a 
tidally  dominant  simulation  is  investigated  using  an  idealized  inlet,  such  as  that  used  by 
Yeeramony  and  Blain  (2001)  and  Kapolnai  et  al.  (1996).  The  idealized  inlet  is  forced  at  the 
offshore  boundary  by  an  M2  tide.  It  is  shown  that  while  the  point-wise  time  series  of 
elevations  are  smoothly  varying  with  minimal  noise,  the  velocity  solutions  are  noisy  and 
have  near  shock-like  characteristics.  These  shock-like  features  are  due  to  the  wetting  and 
drying  of  rows  of  elements  at  a  time,  which  result  in  relatively  large  additions/subtractions  of 
water  mass  to  the  computational  domain.  Even  for  realistic  scenarios,  mesh  nodes  can  be 
aligned  along  bathymetry  contours  or  equivalently  along  the  land/sea  interface.  So  here  we 
examine  how  the  free  parameters  for  wetting  and  drying  impact  the  shock-like  features  that 
can  develop  in  the  velocity  solution. 

For  the  idealized  inlet  shown  in  figure  2.6  the  bathymetry  varies  linearly  from  13 
meters  at  y  =  0.0km  to  approximately  -1.06  meters  above  sea  level  at  y  =  33.0km  and  is 
independent  of  x.  Two  meshes  are  considered  and  shown  in  figure  2.7  along  with  the  location 
of  five  recording  stations  that  transverse  the  land/water  interface.  The  first  mesh  is  highly 
structured,  in  which  regular  divisions  of  equilateral  triangular  elements  are  arranged  in 
regular  rows  and  columns.  These  elements  have  edge  lengths  that  range  from  1000m  at  the 
open  boundary,  to  62.5 m  near  the  opening  of  the  inlet.  For  this  mesh,  the  arrangements  of 
nodes  in  rows  that  follow  bathymetry  contours  results  in  entire  rows  of  elements 
wetting/drying  simultaneously.  To  verify  that  the  resulting  simultaneous 
additions/subtractions  of  large  amounts  of  mass  to  the  simulation  are  the  cause  of  the  shock¬ 
like  features  in  the  velocity  solution,  a  second  unstructured  mesh  is  designed.  For  the  second 
unstructured  mesh,  the  nodes  are  not  arranged  in  orderly  rows  and  columns.  Therefore,  it  is 
not  as  susceptible  to  entire  rows  of  elements  wetting/drying  simultaneously.  This  mesh  has 
roughly  the  same  size  elements  as  the  regular  mesh.  Station  number  1  is  located  on  land, 
station  2  is  at  the  zero  bathymetry  contour,  and  stations  3,  4,  and  5  are  located  in  water  of 
increasing  depth,  respectively.  The  run  is  forced  with  an  M2  tide  of  constant  magnitude 
0.3048m  along  the  open  boundary  for  8  days.  A  4-day  effective  ramp  was  applied  to  the 
forcing  and  a  4  second  time  step  was  required  as  defined  by  the  CFL  condition.  No  other 
forcing  functions  were  applied.  The  hybrid  nonlinear  bottom  friction  option  was  again 
applied  for  this  simulation. 
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Figure  2.6.  The  bathymetry  contours  in  meters  for  the  idealized  inlet  case.  Positive  values  are 
land  values  while  negative  values  are  water. 
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Figure  2.7.  A  close  up  of  (a)  the  structured  mesh  and  (b)  the  unstructured  mesh  showing 
point  locations  (1-5)  for  the  idealized  inlet  case.  The  blue  line  is  the  zero  meter  contour  of 
bathymetry. 

Figure  2.8  shows  the  time  series  of  total  water  depth,  H,  in  meters  and  the  velocity 
components,  u  and  v  in  m/s  at  all  five  stations.  The  first  observation  and  one  most  frequently 
cited  in  others'  works,  i.e.,  Luettich  and  Westerink  (1999),  Dietrich  et  al.  (2004,  2006},  is 
that  the  total  water  depth  varies  smoothly  in  a  cyclic  pattern  following  the  tidal  forcing.  The 
velocity  solutions  are  also  cyclic,  but  they  are  locally  very  noisy.  Our  analysis  will  be 
confined  to  the  dominant  component  of  the  velocity,  the  v-component.  In  figure  2.9,  the  v- 
component  of  velocity  is  examined  over  a  6  hour  period  to  see  details  of  the  (a)  wetting 
phase  and  (b)  the  drying  phase.  During  the  wetting  phase,  stations  3,  2,  and  1  all  become  wet 


18 


in  succession.  Notice  that  just  before  each  station  turns  from  dry  to  wet,  the  v-component  of 
velocity  at  the  stations  behind  it  (in  deeper  water),  begins  to  level  off  or  decrease.  Once  the 
forward  station  wets,  the  deeper  stations  increase  dramatically  in  value  in  a  near  shock-like 
manner.  Stations  in  water  deeper  than  station  5  (not  shown),  and  the  dip  and  surge  pattern 
continues  to  be  exhibited  for  several  rows  of  elements  behind  the  wetting  station.  A  similar 
dip  and  surge  pattern  is  exhibited  during  the  drying  phase  which  is  shown  in  (b). 


Figure  2.8.  Time  series  of  total  water  depth,  H,  horizontal  velocities,  u,v,  at  five  stations. 
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Figure  2.9.  A  6  hour  time  series  of  the  v-component  of  velocity,  at  five  stations  during  (a)  a 
wetting  phase  and  (b)  a  drying  phase.  The  colored  dashed  lines  indicate  when  the  stations 
either  wetted  or  dried. 

The  effects  of  varied  H0  values  on  the  dip  and  surge  characteristics  in  v  and  in  ^ 

at  station  5,  are  shown  in  figure  2.10.  Other  stations  have  similar  results.  The  results  for  the 
unstructured  mesh  are  also  included  for  comparison  purposes.  For  increasing  values  of  H0 , 
the  magnitude  of  the  dip  and  subsequent  surge  becomes  more  severe.  In  comparison,  notice 
how  the  unstructured  mesh  leads  to  a  much  smoother  velocity  solution.  This  is  a  result  of 
fewer  elements  becoming  wet/dry  at  the  same  time,  in  contrast  to  the  row  by  row 
wetting/drying  in  the  uniform  mesh  case.  The  influence  of  the  parameter  Vmin  is  negligible 

on  the  dip  and  surge  feature  observed  in  the  velocity  solution  for  a  fixed  H0  =0.1  m,  as 
shown  in  figure  2.11.  In  fact  the  elevation  curves  for  the  case  of  V  .  =  0.005  m/s  and 

V  .  =0.01  m/s  lay  on  top  of  one  another. 

min  J  r 
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Station  5 


Figure  2.10.  A  6  hour  time  series  at  station  5,  of  the  v-component  of  velocity  and  of  the 
discrete  time  rate  of  change  of  v,  for  different  values  of  H0  over  the  structured  mesh.  The 
velocity  response  over  the  unstructured  mesh  for  H0  =0.01  is  also  included.  The  colored 
dashed  lines  indicate  when  the  station  either  wetted  or  dried  for  the  case  of  H0  =  0.01 . 

The  dip  and  surge  pattern  observed  in  the  v-component  of  velocity,  is  a  result  of  what 
amounts  to  a  temporary  no-normal  flow  boundary  condition  imposed  by  the  WAD  algorithm 
along  the  wet/dry  front.  Along  the  wet/dry  front,  an  artificial  barrier  is  in  place  that  holds  the 
water  back  until  it  reaches  a  critical  height  dictated  by  H0  and  Vmin  .  While  the  barrier  is  in 

place,  water  levels  are  increasing  behind  it  during  a  wetting  phase  and  the  flow  pattern  is 
diverted  due  to  the  no-normal  flow  condition.  Hence,  we  see  a  leveling  off  and/or  decrease 
in  the  v-component  of  velocity  (in  this  case),  until  the  water  level  reaches  the  critical 
threshold  and  the  barrier  is  removed,  there-by  allowing  water  to  surge  forward.  With  the 
regular  mesh,  we  have  a  series  of  temporary  barriers  arranged  in  rows  that  are  encountered 
one  after  the  other  as  inundation  occurs,  but  with  the  unstructured  mesh  this  is  not  the  case. 
During  the  drying  process,  water  levels  decrease  gradually  until  the  cutoff  value  of  Habsm  is 

reached,  after  which  time  the  water  mass  is  removed  from  the  computation.  This  sudden 
removal  of  water  mass  causes  the  shock-like  condition  seen  in  the  drying  phases  in  the  v- 
component  of  velocity.  The  height  of  the  artificial  barriers  is  controlled  principally  by  H0 
and  that  is  why  we  see  the  most  influence  on  the  dip  and  surge  pattern  from  varying  H0 . 
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Figure  2.11.  A  6  hour  time  series  of  the  v-component  of  velocity,  at  station  5  for  a  wetting 

phase  with  varied  values  of  V  .  .  Note  that  the  results  for  cases  V  •  =0.01  and 
^  min  min 

Vmin  =  0.005  are  virtually  indistinguishable.  The  colored  dashed  lines  indicate  when  the 

stations  wetted  for  the  case  of  V  •  =  0.01 . 

min 


The  variability  of  ^ ,  the  acceleration,  demonstrates  the  influence  of  the 


wetting/drying  of  nearby  elements  on  the  velocity  solution.  Of  particular  importance  is  how 
noisy  the  values  of  ^  are  immediately  following  a  wetting/drying  of  a  row  of  elements. 


Not  only  are  the  velocity  solutions  experiencing  sudden  changes  but  so  too  are  the  values 
of£ .  These  sudden  changes  can  lead  to  instabilities  during  the  solution  process.  The 


authors'  own  experiences  confirm  that  instabilities  in  the  WAD  are  almost  always  first 
observed  in  the  velocity  solution  step. 


d.  Gently  Sloping  Beach 

While  the  wetting  and  drying  algorithm  performs  well  in  tidal  regions  that  have 
simple  bathymetric  slopes,  there  are  limitations  when  the  bathymetric  slope  flattens,  as  in  the 
case  of  a  tidal  or  mud  flat.  The  Gently  sloping  Beach  test  problem  demonstrates  that  the 
algorithm  allows  wetting  and  drying  fronts  to  propagate  in  non-symmetric  ways,  even  for 
symmetric  applications  such  as  this  one.  Additionally,  during  the  drying  phase,  small 
pockets  of  elements  can  remain  wet  even  when  negative  pressure  gradients  exist.  The  mesh 
for  the  gently  sloping  beach  is  designed  to  highly  resolve  the  wetting  and  drying  region  and 
to  smoothly  transition  to  decreasing  element  sizes  with  decreasing  water  depth,  see  figure 
2.12a.  The  spatial  scales  of  the  problem  are  proportional  to  those  of  a  realistic  beach  under 
tidal  forcing.  The  modeled  domain  is  a  square  box  of  dimension  4,960  meters  with  a 
bathymetry  profile  that  gently  slopes  upward  to  a  plateau  beach.  The  water  depth  ranges 
from  9.32 m  at  the  open  ocean  boundary  to  -0.1m  above  sea  level  at  the  plateau  beach.  The 
bathymetry  profile  in  the  cross-shore  (along  the  y-axis)  is  as  follows, 
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2.0-  0.0060  - 1 200), if  y  <  1 200m; 


z 


(y)=\ 


co  +ciy+c2y2  +c3j3,if  y  <  2500 m; 


— 0.1,  if  jv  >  2500 m 


where  the  cubic  polynomial  coefficients  c,  are  determined  by  solving  a  4  x  4  system  of  linear 

dz 

equations  satisfying  z(l  220m)  =  2.0 m ,  z(1780m)  =  -0. 1  m  and  the  slopes  (z'  =  —  )  satisfying, 

dy 

z'(1220m)  =  -0.006  and  z'(1780m)  =  0 m  .  The  mesh  consists  of  92,176  computational  nodes 
and  183,836  triangular  elements.  The  mesh  spacing  is  graded  from  the  open  ocean  boundary 
where  nodal  spacing  is  80m  to  the  water/land  transition  zone  (bathymetry  +/-)  where  nodal 
spacing  is  only  5m,  see  figure  2.12a. 

The  model  simulation  is  for  14.4  hours  with  a  6  hour  ramp  period  applied  to  the 
forcing.  A  0.1  second  time  step  is  required  by  the  CFL  condition.  The  simulation  is  forced 
for  9  hours  with  an  elevation  specified  as  a  sinusoidal  wave  with  amplitude  of  0.135  meters 
and  a  period  of  one  hour.  After  9  hours  the  tidal  forcing  function  is  shut  off  entirely  and  no 
other  forcing  is  applied.  The  period  without  forcing  allows  for  investigation  of  the  drying 
phase  of  the  WAD  algorithm.  A  value  of  H0  =  0.01  meters  for  the  minimum  water  depth  and 
a  value  of  Vmin  =  0.01  m/s  for  the  minimum  wetting  velocity  are  selected. 

During  the  tidally  forced  portion  of  the  run,  the  wetting  and  drying  front  moves  in  a 
predictable  fashion.  Water  moves  up  the  slope  and  onto  the  beach  plateau  during  inundation. 
This  is  subsequently  followed  by  a  retreat  of  the  wetting  front  off  the  beach  and  down  the 
slope.  An  unexpected  result  is  the  non-symmetric  propagation  of  the  wetting  front,  see  figure 
2.12b.  This  same  asymmetry  is  observed  during  the  drying  phase.  A  matter  of  more  concern 
is  the  lack  of  drainage  during  quiescent  (non-forced)  periods,  see  figures  2.12c  and  2.12d. 
The  pockets  of  wet  elements  left  on  the  beach  remain  active  in  the  simulation  and  can 
contribute  to  instabilities.  In  a  wind-driven  application,  typical  for  inundation  concerns,  such 
instabilities  can  easily  manifest  when  pockets  of  water  move  around  in  response  to  rising 
elevations  caused  by  the  wind  forcing.  In  fact,  since  the  WAD  algorithm  is  not  mass 
conserving,  it  would  be  possible  to  cause  pockets  of  water  to  either  increase  in  area,  stay  at  a 
constant  depth  or  even  increase  in  depth,  even  though  these  pockets  are  cut  off  from  the  main 
body  of  water. 

As  we  mentioned  earlier,  one  of  the  past  upgrades  to  the  WAD  algorithm  was  an 
empirical  check  that  prevented  downhill  flow  from  originating  from  barely  wet  nodes,  except 
in  the  case  of  overtopping  a  weir  (or  barrier  node).  The  concept  was  conceived  for  situations 
when  water  levels  crept  up  a  hill  and  barely  wet  the  crest  of  the  hill  which  would  in  turn 
cause  downward  flow  behind  the  hill.  However,  the  same  scenario  is  occurring  for  this 
gently  sloping  plateau  beach  case  when  we  are  trying  to  drain  water  back  off  the  plateau.  In 
this  case,  the  empirical  check  is  turning  off  elements  (drying  them),  before  all  the  water  is 
drained  and  leaving  pockets  of  wet  elements,  see  for  example  figure  2.12c  between 
y  =  1670m  and  y  =  1725m. 
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e.  Summary 


A  detailed  sensitivity  analysis  of  the  surge  and  current  response  to  the  user  controlled 
parameters  of  the  wetting  and  drying  algorithm,  H0  and  Vmm ,  used  in  the  continuous 

Galerkin  formulation  of  ADCIRC  has  been  presented  Three  case  studies  support  the 
analyses  and  bring  to  light  some  computational  issues  observed  with  the  WAD  scheme. 
Issues  include  a)  noisy  velocity  solutions  caused  by  the  sudden  wetting  or  drying  of  large 
numbers  of  elements  simultaneously,  b)  difficulty  in  drying  elements  over  very  small 
bathymetric  gradients,  a  problem  for  tidal  flats,  and  c)  instabilities  that  result  from 
incomplete  drainage  leaving  isolated  pockets  of  wet  elements.  This  last  issue  is  very 
problematic  during  wind-forced  simulations  in  that  the  wind  acts  on  these  isolated  elements 
of  water.  The  parameter  Ha  is  most  influential,  affecting  the  global  mass  balance  due  to  the 
non-conservative  nature  of  the  WAD  algorithm.  The  parameter  Vmin  has  substantially  less 

effect  on  the  WAD's  performance.  Depending  on  the  purpose  of  the  simulation,  a  small 
value  of  H0  does  not  always  produce  the  best  solution.  However,  larger  Ha  values  can 
produce  shock-like  characteristics  in  the  velocity  solutions.  Overall,  for  most  cases  tested  the 
WAD  algorithm  produces  reasonable  results  but  care  should  be  taken  in  setting  the  value  of 
H0  .  One  side  note,  using  external  boundary  types,  IBTYPE  =  30,  40,  and  41  with  boundary 
nodes  that  are  allowed  to  wet  and  dry,  can  quickly  result  in  an  unstable  velocity  solution 
under  cyclic  forcing  conditions.  This  is  due  in  part  to  the  fact  that  no  constraints  are  being 
placed  on  the  velocity  solution  for  these  boundary  types  and,  as  such,  the  velocities  can 
become  unrealistically  large  during  the  ebb  stage  of  the  cycle. 
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Figure  2.12.  For  the  gently  sloping  beach,  (a)  shows  a  close  up  of  the  mesh 
topography/bathymetry  with  nodes  indicated  by  blue  circles,  the  magenta  diamond  indicates 
the  last  node  below  MSL  and  the  green  diamond  indicates  the  first  node  on  the  plateau  beach, 
-0.1m  above  sea  level.  In  (b)-(c)  a  series  of  snap  shots  of  elevation  (dry  nodes  indicated  with 
white  ’.’)  for  (b)  a  wetting  phase,  (c)  a  drying  phase  and  (d)  the  final  elevation. 


4.  Validation  Test  Results  -  Hurricane  Katrina  Storm  Surge 

Throughout  modem  history,  amphibious  assaults  and  landings  have  been  a  mainstay 
of  U.S.  Navy  operations.  The  vulnerability  of  these  landing  craft  to  capsizing,  swamping, 
stranding,  and  filling  with  sand  and  water  was  clearly  realized  following  a  post-World  War  II 
review  of  amphibious  operations.  Many  amphibious  landing  problems  and  casualties  during 
World  War  II  could  be  attributed  to  the  waves,  currents  and  water  levels  of  the  local 
environment.  Following  the  major  invasion  of  Incheon  Harbor  in  the  Republic  of  Korea,  a  U. 

5.  Navy  Tank  Landing  Ship  was  stranded  during  low  tide  near  the  Tidal  Basin  on  Incheon’s 
waterfront,  20  September  1950. 
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More  than  fifty  years  later,  the  Navy  still  finds  inundated  environments  challenging 
for  operations.  Since  the  declaration  of  the  Global  War  on  Terrorism  following  the  events  of 
September  11,  2001,  military  operations  are  increasingly  focused  on  special  operations  that 
take  place  in  coastal  environments  such  as  estuaries,  shallow  waterways,  and  inland  rivers. 
The  occurrence  of  inundation  in  these  operational  theatres  is  typically  caused  by  extreme 
tidal  ranges,  rainfall-induced  flooding  events,  and/or  wind-generated  setup  that  directly  affect 
the  insertion  and  movement  of  Naval  Special  Warfare  (NSW)  forces.  NSW  forces  routinely 
insert  themselves  into  river  and  estuarine  environments  where  the  only  known  information 
may  be  an  outdated,  perhaps  30-year  old,  topographic  map.  Areas  that  are  subject  to 
inundation  processes  are  often  located  at  the  cusp  of  the  land-sea  interface  where  algorithms 
for  processing  satellite  imagery  break  down  or  are  sub-optimal. 

Inundation  from  storm  surge  is  also  a  concern  for  stateside  Navy  installations.  The 
two  major  homeports  for  the  U.  S.  Navy’s  east  coast  fleet  are  at  Norfolk,  Virginia  and 
Mayport,  Florida,  both  vulnerable  to  landfalling  Atlantic  hurricanes.  A  decision  to  relocate 
the  Norfolk  harbor  fleet,  for  example,  could  cost  $5  million  and  would  need  to  take  place 
three  days  in  advance  of  a  predicted  landfall  in  order  to  recall  personnel  and  ready  ships  in 
maintenance  or  overhaul  for  evacuation.  Most  recently,  the  Navy  base  at  Pascagoula, 
Mississippi  located  on  the  Gulf  of  Mexico  was  directly  impacted  by  the  landfall  of  Hurricane 
Katrina,  29  August  2005. 

As  we  now  know,  the  Naval  Station  Pascagoula  was  not  alone  in  registering  effects 
from  Hurricane  Katrina.  The  devastation  to  Gulf  Coast  communities  on  29  August,  2005 
from  Hurricane  Katrina  far  exceeded  all  previously  recorded  storm  events.  From  the  extent  of 
storm  damage  to  the  coastal  states  of  Louisiana,  Mississippi  and  Alabama,  Katrina  has  been 
categorized  as  the  most  destructive  and  costliest  natural  disaster  in  the  history  of  the  United 
States.  According  to  the  National  Oceanic  and  Atmospheric  Administration  (NOAA),  the 
storm  surge  along  the  Mississippi  coast  was  the  highest  storm  surge  ever  recorded  in  the 
United  States. 

The  storm  surge  and  inundation  from  Hurricane  Katrina  that  devastated  Mississippi 
Gulf  Coast  communities  on  29  August,  2005  presents  a  unique  opportunity  to  evaluate  the 
capabilities  at  NRL  and  within  the  Navy  to  predict  storm  surge  and  inland  inundation. 
Reconstruction  of  the  storm  surge  and  inundation  events  precipitated  by  Hurricane  Katrina 
provide  an  invaluable  opportunity  to  evaluate  the  Navy’s  capability  to  predict  coastal  surge 
and  inundation  and  to  direct  future  developments  that  enhance  such  a  capability.  A  highly 
realistic  simulation  of  Katrina’s  storm  surge  and  inland  inundation  is  developed  using  the 
2DDI  version  of  the  ADvanced  CIRCulation  (ADCIRC)  model.  Unprecedented  observations 
of  the  currents  as  well  as  high  water  marks,  extent  of  inland  inundation  and  water  levels  are 
available  to  assess  model’s  performance. 

a.  Model  Configuration 

The  initial  requirement  for  reconstruction  of  hurricane  Katrina’s  storm  surge  along 
the  Mississippi  Gulf  Coast  was  a  computational  mesh  that  extended  from  the  shoreline  to 
inland  locations.  The  importance  of  a  quality  mesh  cannot  be  understated  in  the  modeling  of 
inundation  events.  To  accurately  represent  the  surge  and  inundation  the  mesh  must  resolve 
fine-scale  changes  in  bottom  slope  and  topography,  details  of  the  coastline,  and  other 
geographic  features  such  as  islands,  inlets  and  channels,  while  simultaneously  preserving 
properties  of  the  triangular  elements  that  promote  model  stability  and  retaining  a 
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computationally  viable  problem  (i.e.  a  timely  solution).  The  transitioned  software,  MeshGUI, 
developed  for  semi-automated  mesh  generation  is  applied  to  construct  an  unstructured  finite 
element  mesh  using  refinement  criteria  based  on  the  specified  bathymetric/topographic 
values  and  constrained  by  the  coastal  boundary  points  that  extend  over  the  land  to  allow 
inundation  processes.  The  bathymetric  and  topographic  data  source  used  in  shallow  coastal 
areas  was  taken  from  Northern  Gulf  Littoral  Initiative  bathymetric  data  at  3  arc  second 
resolution.  To  eliminate  the  discontinuities  between  data  blocks  within  the  NGLI  data  set, 
particularly  in  the  Pearl  River  Basin,  bilinear  interpolation  was  applied  across  such 
boundaries  to  smooth  the  bathymetry  and  topography  values.  Smoothly  varying  bathymetry 
and  topography  from  element  to  element  is  very  important  for  the  stability  of  ADCIRC’s 
wetting  and  drying  algorithm.  The  final  unstructured  finite  element  mesh  designed  to  best 
capture  Katrina’s  storm  surge  and  inland  inundation  consists  of  489,071  nodes  and  956,869 
triangular  elements  (figure  4.1).  Resolution  ranges  from  63 m  to  nearly  6km  (figure  4.2)  with 
225 m  to  500m  resolution  in  most  shallow  coastal  and  inland  areas. 

The  mesh  centers  on  the  northern  Gulf  Coast  region  encompassing  inland  areas,  but 
also  includes  the  entire  Gulf  of  Mexico  and  extends  out  into  the  western  North  Atlantic 
Ocean.  Such  an  expansive  domain  allows  the  surge  to  naturally  build  up  within  the  modeled 
region  as  the  hurricane  moves  from  the  deep  ocean  into  coastal  waters  (Blain  et  al.,  1994). 
Furthermore,  ocean  boundaries  in  deep  water  are  subject  to  minimal  surge  and  inverted 
barometer  effects  and  can  appropriately  accept  tidal  forcing  from  a  global  tide  model.  These 
boundaries  also  are  removed  from  the  coastal  area  of  interest. 

The  landward  boundary  includes  inland  areas  concentrated  on  the  Mississippi  coast 
extending  eastward  just  past  Mobile  Bay,  AL  and  westward  including  the  entire  Pearl  River 
basin  and  areas  inland  to  Interstate  10  on  the  north  shore  of  Lake  Ponchartrain,  LA.  As 
important  conduits  of  storm  surge,  the  inclusion  of  river  basins  from  the  coast  to  inland, 
upstream  locations  is  important.  The  inundating  surge  did  not  actually  penetrate  to  the  land 
boundaries  rendering  the  no  flow  condition  (IBTYPE  =  0)  a  good  choice.  If  surge  had 
propagated  to  the  inland  boundary,  a  radiative  boundary  condition  would  be  more 
appropriate  (IBTYPE  =  30)  to  prevent  the  surge  from  reflecting  back  into  the  model  domain. 
The  targeted  spatial  resolution  of  the  mesh  near  the  coast  and  inland  is  225 m  and  represents  a 
balance  between  the  desire  for  fine-scale  resolution,  the  need  for  stability  of  the  inundation 
algorithm,  and  accounts  for  computational  constraints  imposed  by  the  necessarily  small  time 
step  integration. 
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Figure  4.1.  The  density  of  the  triangular  elements  in  the  northeast  Gulf  of  Mexico  contained 
in  the  computational  mesh  for  the  entire  model  domain  (inset)  overlays  a  true  color  satellite 
image  that  distinguishes  land  from  water. 
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Figure  4.2  The  resolution  in  meters  of  the  computational  mesh  created  for  the  northeast  Gulf 
of  Mexico  storm  surge  and  inundation  computations. 

To  drive  the  surge  model,  the  best  available  wind  forcing  was  produced  by  NOAA’s 
Hurricane  Research  Division  (HRD)  at  the  Atlantic  Oceanographic  and  Meteorological 
Laboratory  (AOML)  through  the  HRD  Real-time  Hurricane  Wind  Analysis  System 
(H*Wind)  project.  The  H*Wind  product  is  an  integrated  tropical  cyclone  observing  system 
within  which  wind  measurements  from  a  variety  of  observation  platforms  are  used  to  develop 
an  objective  analysis  of  the  distribution  of  wind  speeds  in  a  hurricane  (Powell  et  al.,  1998). 
The  wind  fields  are  typically  constructed  from  a  real-time  analysis  of  flight-level 
reconnaissance  data,  satellite  observations,  pressure-wind  relationships  and  available  surface 
data.  We  interpolate  the  3-hourly  H*Winds  using  an  approach  that  follows  the  storm  center 
in  time  to  preserve  the  integrity  of  the  storm  as  it  moves  in  time  and  we  further  downscale  the 
winds  fields  to  fifteen-minute  intervals.  The  time  interpolated  wind  fields  are  spatially 
interpolated  to  the  computational  mesh  and  then  converted  to  wind  stress.  The  wind  drag  at 
the  sea  surface  is  simply  specified  as  a  constant  and  with  no  distinction  between  winds  over 
land  or  water  or  the  directional  history  of  the  wind.  More  recent  versions  of  the  ADCIRC 
model  now  accommodate  spatially  varying  surface  roughness  values  that  can  be  computed  to 
reflect  land  use,  canopy,  and  vegetative  cover  leading  to  more  realistic  representations  of  the 
wind  drag  (Westerink  et  al.,  2008). 

In  addition  to  surface  winds,  tidal  forces  are  applied  including  those  that  act  on  the 
modeled  body  of  water  (tidal  potential)  and  those  caused  by  tides  entering  the  domain  at  the 
open  ocean  boundary.  At  the  deep  ocean  boundary,  tidal  forcing  is  applied  at  frequencies  of 
the  daily  (Ki,  Oi)  and  twice  daily  (M2,  S2,  and  N2)  tides  obtained  from  the  Grenoble  global 
tidal  model  (FES99,  Lefevre  et  al,  2002).  The  tidal  potential  is  applied  on  the  interior  of  the 
domain  for  the  same  constituents. 
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The  ADCIRC  model  fort.  1 5  input  file  parameter  specifications  for  the  2D  wetting  and 
drying  application  with  tide  and  wind  forcing  are  configured  using  the  Makefl5  GUI.  The 
default  settings  in  the  Makef  1 5  GUI  were  changed  in  the  following  manner: 

■  Activate  the  non- fatal  error  override,  NFOVER  =  1 

■  A  hybrid  nonlinear  bottom  friction  formulation  is  selected,  NOLIBF  =  2 

■  Activate  tidal  potential  forcing,  NTIP  =  1 

■  Wind  velocity  and  surface  pressure  forcing  are  selected,  NWS  =  4 

■  Time  step  is  reduced  to  1  sec,  DTDP  =  1.0 

■  The  reference  time  is  set  to  the  time  after  ramping,  REFTIM  =  14.0 

■  The  meteorological  forcing  time  increment  is  15  minutes  (900  sec), 
WTIMINC  =  900 

■  The  length  of  simulation  is  set  to  18.41666667  days,  RNDAY  =  18.41666667 

■  The  ramp  period  is  set  to  a  duration  of  15  days,  DRAMP  =15.0 

■  The  minimum  depth,  HO,  is  set  to  0.01  and  the  minimum  velocity  for  wetting 
is  0.01  (NODEDRYMIN  and  N ODE WETMIN  are  not  used),  HO  =  0.01, 
VELMIN  =  0.01. 

■  The  central  projection  points  are  specified  for  the  model  domain, 
SLAM0,SFEA0  =  -88.5,  29.0,  respectively. 

■  Parameters  for  the  hybrid  nonlinear  bottom  friction  coefficient  are  specified, 
CF  =  0.003  (minimum  friction  coefficient),  HBREAK  =  2.0  (break  depth), 
FTHETA  =  10.0,  and  FGAMMA  =  1.333  (see  equation  (6),  Section  lb) 

■  Assign  5  tidal  potential  constituent,  NTIF  =  5 

■  Tidal  potential  constituents,  TIPOTAG  =  Kl,  01,  M2,  S2,  N2 

■  Input  the  date  at  the  start  of  the  simulation  to  compute  the  nodal  factors 
(August  26,  2005) 

■  Assign  6  periodic  boundary  forcing  frequencies,  NBFR  =  6 

■  Tidal  potential  constituents,  BOUNTAG  =  STEADY,  Kl,  01,  M2,  S2,  N2 

■  Input  the  date  at  the  start  of  the  simulation  to  compute  the  nodal  factors 
(August  26,  2005) 

■  Select  global  elevation  output,  NOUTGE  =  -1  for  just  under  3  days 
(TOUTSGE  =  15.97916667,  TOUFGE  =  18.416667.0  model  days)  every  half- 
hour  (NSPOOLE  =  1800) 

■  Select  global  velocity  output,  NOUTGY  =  -1  for  just  under  3  days 
(TOUTSGV  =  15.97916667,  TOUFGV  =  18.416667.0  model  days)  every 
half-hour  (NSPOOLV  =  1800) 

■  Output  a  hotstart  every  6  hours  (21600  time  steps),  NHSTAR  =  1,  NHSINC  = 
21600  time  steps 

The  hindcast  simulation  of  Hurricane  Katrina  storm  surge  began  27  August  2005 
(0000  UTC)  following  a  ramp-up  period  of  15  days  during  which  all  forcings  were  gradually 
applied  until  full  strength  was  reached  at  the  end  of  the  ramp-up  phase.  By  this  time  in  the 
simulation,  Hurricane  Katrina  had  crossed  the  state  of  Florida  and  had  entered  the  warm 
waters  of  the  Gulf  of  Mexico  (figure  4.3).  Katrina  was  past  its  peak  intensity  by  the  first 
landfall  near  Buras,  Louisiana  at  6:10  am  CDT  (1110  UTC)  29  August  2005  and  a  second 
landfall  near  the  Louisiana/Mississippi  border  occurred  about  9:45am  CDT  (1445  UTC)  29 
August  2005.  The  model  hindcast  of  surge  and  inundation  ended  at  5:00am  CDT  (1000 
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UTC)  on  30  August  2005  which  coincided  with  the  last  available  H*wind  product  contained 
within  the  mesh.  At  every  1-sec  time  integration  of  the  model,  the  water  levels  and  depth- 
integrated  currents  are  computed  by  the  ADCIRC  model  at  all  points  in  the  model  domain. 


Figure  4.3:  The  track  of  Hurricane  Katrina  at  6  hour  intervals  once  the  storm  enters  the  Gulf 
of  Mexico  off  the  west  Florida  coast  at  6Z  26  August  2005  with  central  pressure  of  987  mb. 
The  final  location  of  the  storm  at  hurricane  strength  occurred  in  southern  Mississippi  at  18Z 
29  August  2005  with  a  central  pressure  of  948  mb. 

b.  Water  Levels 

The  surface  winds  from  Hurricane  Katrina  at  10:00am  CDT  (1500  UTC)  on  29 
August  2005  (figure  4.4)  and  the  resulting  storm  surge  computed  by  the  ADCIRC  model 
(figure  4.5)  for  the  same  time  and  date  are  presented.  Surge  heights  well  over  20  feet  at  the 
coastline  on  the  right  side  of  storm  reflect  not  only  the  strength  of  the  storm  winds  at  the  time 
of  landfall  but  also  the  build-up  of  surge  that  occurred  prior  to  landfall.  Even  higher  water 
levels  are  shown  inland  (near  Waveland  Mississippi  and  west  of  Biloxi,  Mississippi)  as  the 
large  radius  of  hurricane  winds  easily  pushed  water  over  the  gently  sloping  coastal  lands. 
Timing  of  the  inundation  indicates  that  areas  to  the  west  of  Waveland,  Mississippi  including 
the  northern  coast  of  Louisiana  (Slidell)  inundated  first  as  hurricane  winds  pushed  water  into 
the  bays  and  up  the  rivers.  Not  until  landfall  did  the  Mississippi  Gulf  coast  west  of  Biloxi 
experience  its  peak  flooding.  Note  that  even  after  landfall  (figure  4.5),  sea  levels  remain 
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elevated  throughout  the  coastal  waters.  For  some  areas  far  inland  particularly  at  the  wetting 
front,  excessive  inland  flooding  (over  30  feet)  is  computed.  Analysis  of  these  hindcast  results 
have  revealed  limitations  in  the  inundation  mechanism  within  the  ADCIRC  model  that 
prevent  rapid  advancement  of  a  wetting  front  and  incomplete  drainage  of  the  flood  water 
following  peak  storm  winds  (evidence  of  this  is  northwest  of  Stennis  Space  Center, 
Mississippi). 

Timing  of  the  wetting  front  is  difficult  to  validate  since  observations  are  often  limited. 
However,  the  modeled  water  heights  are  compared  to  recorded  elevations  at  three  observing 
stations  that  survived  the  storm,  Pilot’s  Station,  SW  Pass,  Louisiana,  Waveland,  Mississippi 
and  Dauphin  Island,  Alabama  (figure  4.6).  In  each  case  the  agreement  between  the  modeled 
and  observed  water  levels  is  quite  reasonable  with  correlation  coefficients  of  0.75  or  higher. 
The  phasing  of  the  tides  and  peak  surge  computed  by  the  model  are  leading  the  observed 
values  by  no  more  than  a  couple  of  hours  and  the  model  underprediction  as  the  storm  nears 
its  landfall  position  is  likely  due  to  the  neglect  of  surface  wave  effects. 
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Figure  4.4:  The  magnitude  (color)  and  direction  (arrows)  of  the  maximum  1 -minute  sustained 
surface  winds  in  knots  for  Hurricane  Katrina  at  10:00am  CDT  (1500  UTC)  on  29  August 
2005.  [Courtesy  of  the  NOAA  Hurricane  Research  Division] 
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Elevation  (ft)  on  29- August-2005  at  10:00:00  AM 
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Figure  4.5:  The  ADCIRC  model  computed  storm  surge  and  inland  inundation  elevation  in 
feet  for  Hurricane  Katrina  at  10:00am  CDT  (1500  UTC)  on  29  August  2005.  The  coastal 
outline  is  shown  in  black.  The  storm  center  is  shown  by  a  circled  X  and  the  location  of 
Stennis  Space  Center,  Mississippi  is  given  by  a  star. 


33 


Station  ID:  8735180  —  Dauphin  Island,  AL 


Figure  4.6:  Time  series  of  model  computed  (blue)  and  observed  (red)  water  elevations  in  feet 
at  three  NOAA  coastal  stations,  Pilots  Station  East,  Sw  Pass,  Louisiana,  Waveland 
Mississippi,  and  Dauphin  Island,  Alabama. 

c.  Inundation 

A  qualitative  assessment  of  the  model-computed  inundation  is  provided  by  comparing 
a  map  of  modeled  high  water  levels  with  Federal  Emergency  Management  Agency  (FEMA) 
inundation  maps  (FEMA,  2005  and  FEMA,  2006)  as  shown  in  figure  4.7.  On  the  FEMA 
maps  for  Louisisana  (figure  4.7b)  and  Mississippi  (figure  4.7c),  zones  of  inundation  are 
highlighted  in  yellow  with  the  storm  track  visible  as  a  red  line  on  the  Louisiana  map.  The 
FEMA  and  ADCIRC  inundation  maps  are  alligned  vertically  within  figure  4.7  to  facilitate  a 
qualitative  intercomparison  between  the  two.  The  ADCIRC  inundation  map  (figure  4.7a)  is 
created  by  recording  the  highest  water  elevation  at  each  model  grid  point  over  the  entire 
simulation  period  of  the  storm.  These  high  water  mark  maps  provide  no  information  on 
timing  of  the  inundation  front  and  the  comparison  to  FEMA  inundation  maps  does  not 
provide  an  assessment  of  the  magnitude  of  the  water  level,  but  rather  this  comparison 
measures  the  model  capability  with  repsect  to  the  spatial  extent  of  inland  inundation. 
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Moving  from  west  to  east  in  figure  4.7,  the  model  represents  well  the  significant 
pathway  for  inland  inundation  in  Lousiana  along  the  Pearl  River  basin.  The  lengthly  inland 
extent  of  inundation  associated  with  the  Pearl  River  basin  indicated  in  the  FEMA  map  is  also 
well  captured  by  the  model  with  notably  less  inundation  to  the  west  of  the  Pearl  River  basin. 
However,  the  model  records  anomalous  inundation  just  to  the  east  of  the  hurricane  Katrina 
landfall  path  in  Mississippi.  The  FEMA  inundation  map  for  Mississippi  (figure  4.7c)  shows 
large  areas  of  inundation  west  of  Bay  St.  Louis  which  is  consistent  with  locations 
immediately  to  the  right  side  of  a  landfalling  hurricane.  Moving  eastward  the  primary 
inundation  pathways  are  the  tributaries  entering  Bay  St.  Louis  and  Biloxi  Bay.  Very  little 
inundation  occurs  along  the  coastline  between  Bay  St.  Louis  and  Biloxi  Bay  which  is 
indicative  of  locally  high  ground.  The  modeled  inundation  also  reflects  this  pattern  of 
inundation  (following  river  tributaries  out  of  bays  and  limited  inland  inundation  along  the 
coast  near  Pass  Christian,  MS).  The  exception  in  the  modeled  inundation  is  the  more 
extensive  inundation  casued  by  the  merging  of  water  from  the  eastern  tributaries  of  Bay  St. 
Louis  with  water  coming  from  the  western  tributaries  of  Biloxi  Bay.  The  remaining  area  of 
significant  inundation  both  recorded  by  the  FEMA  map  and  computed  by  the  model  is  in  and 
around  Pascagoula  Bay  and  inland  up  the  Pascagoula  River.  Overall  the  modeled  patterns  of 
inundation  reflect  those  shown  in  the  FEMA  Inundation  Maps.  In  the  model,  it  is  likely  that 
the  wetting  front  moved  inland  too  quickly  covering  more  area  in  the  regions  of  highest 
water  levels  and  then  steeper  topography  prevented  additional  wetting  and  caused  water  to 
pile  up  artifically  in  several  locations  as  inidated  by  the  red  in  the  figure  4.7a. 

An  evaluation  of  the  modeled  high  water  magnitudes  of  the  surge  and  inundation  is 
accomplished  by  comparing  computed  high  water  mark  values  to  high  water  marks  measured 
by  the  United  Stated  Geological  Survey  shortly  after  the  storm.  At  each  location  in  the  mesh, 
the  highest  water  level  from  the  model  (evaluated  on  10-minute  intervals)  is  recorded  and 
shown  in  Figure  4.8a.  Of  458  high  water  mark  stations,  315  were  wetted  in  the  model.  Red 
dots  on  the  map  in  figure  4.8a  indicate  143  locations  that  did  not  experience  inundation 
during  the  hindcast  simulation  of  Hurricane  Katrina.  It  is  likely  that  a  number  of  factors 
contribute  to  this  type  of  error,  i.e.,  erroneous  values  for  local  water  depth  and  land  height, 
not  accounting  for  the  decreased  wind  drag  over  water,  or  limitations  in  the  inunndation 
mechanism  of  the  model  as  previously  documented  in  Section  3.  Despite  the  non- wetting  of 
certain  locations,  the  model  computed  water  elevations  at  the  remaining  3 1 5  high  water  mark 
locations  had  an  average  error  of  only  1.2  feet  (figure  4.8b).  Stations  with  the  largest  errors 
underpredict  water  levels  and  are  found  near  those  same  locations  that  remained  erroneously 
dry.  Note  that  dry  area  locations  coincide  with  the  limited  inland  inundation  shown  in  the 
FEMA  map  of  figure  4.7c.  A  spatial  resolution  in  the  model  of  225 m  may  have  been  too 
coarse  to  accuractly  capture  the  wetting  front  in  this  region.  Overall,  the  exhibited  model  skill 
is  extraordinary,  given  that  the  hindcast  only  used  readily  available  information  on  water 
depth,  land  height,  and  wind  strength.  The  level  of  detail  reflects  typical  conditions  for  Navy 
operations  in  non-US  waters. 
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Computed  High  Water  Level  ADCIRC  Model 


FEMA  Inundation  Maps 


Figure  4.7:  (a)  A  map  of  ADCIRC  model-computed  high  water  marks  for  Hurricane  Katrina; 
(b)  the  FEMA  Louisiana  Hurricane  Katrina  Surge  Inundation  Map  and  c)  the  FEMA 
Mississippi  Hurricane  Katrina  Surge  Inundation  Map. 
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USGS  High  Water  Mark  Comparisons 
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Figure  4.8:  (a)  A  map  of  inundated  (blue)  or  dry  (red)  U.S.G.S  high  water  mark  locations  as 
computed  by  the  model;  (b)  Comparisons  of  the  modeled  (blue)  and  measured  (red)  high 
water  elevations  in  feet  at  315  USGS  stations. 


d.  Currents 

The  University  of  Southern  Mississippi  had  deployed  Buoy  USM3m001  as  an  initial 
observing  element  of  the  Central  Gulf  of  Mexico  Ocean  Observing  System  on  December  13, 
2004  at  the  location  30°  02’  32.710”N  and  88°  38’  50.235”W  (figure  4.9a).  Data  were 
telemetered  every  3  hrs  via  Globalstar  until  it  was  recovered  after  hurricane  Katrina  in 
September  of  2005.  The  eye  of  hurricane  Katrina  passed  49 nm  to  the  west  of  the  3m  discus 
buoy,  in  20  water  depth,  operated  by  the  Central  Gulf  of  Mexico  Ocean  Observing  System. 
The  buoy  moved  2 km  to  the  northwest  during  the  storm  surge  and  then  13 km  to  the  southwest 
as  the  surge  retreated  (figure  4.9b).  The  buoy  begins  to  drag  its  mooring  when  the  significant 
wave  height  (SWH)  exceeds  6m  (A  in  figure  4.10).  Buoy  latitude  reaches  a  maximum  when 
the  eye  of  Katrina  shares  the  same  latitude  as  the  nominal  buoy  position  (B  in  figure  4.10). 
Southward  movement  of  the  buoy  stops  when  the  significant  wave  height  falls  below  4m  (C 
in  figure  4.10).  Waves  built  to  nearly  half  the  mean  water  depth  and  wave  height  rapidly 
decayed  5-6  hours  after  the  hurricane  moved  over  land. 
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Elevation  (ft)  on  29- August- 2005  at  10:00:00  AM 
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Figure  4.9:  (a)  The  location  of  Buoy  USM3m001  shown  as  a  blue  dot  and  (b)  the 
northwestward  movement  of  the  buoy  during  hurricane  Katrina,  August  2005. 
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Figure  4.10:  Time  series  of  the  latitude  of  Katrina’s  eye,  the  latitude  of  Buoy  USM3m001 
and  the  buoy  measured  significant  wave  heights.  At  A  the  buoy  begins  to  drag  the  mooring 
when  SWH  >  6  m;  at  B,  buoy  latitude  reaches  a  maximum  when  Katrina’s  eye  is  at  the  same 
latitude  as  nominal  buoy  position;  at  C,  the  buoy  stops  moving  southward  when  SWH  <  4  m. 

Figure  4.11  compares  the  measured  buoy  wind  stress  to  the  computed  stress  as 
applied  to  the  model  at  the  location  of  the  buoy.  The  wind  stress  applied  to  the  model  over¬ 
predicts  the  east-west  component  of  the  wind  with  the  most  significant  over-prediction 
occurring  just  prior  to  landfall  (figure  4.11a).  The  northward  component  of  the  wind  stress 
applied  to  the  model  is  slightly  under-predicted  at  the  peak  but  remains  artificially  elevated 
after  the  peak  winds  have  passed  (figure  4.1  lb). 

Currents  show  an  asymmetrical  storm  surge  event  with  stronger  currents  during  the 
fall  of  the  surge  (figure  4.11).  Comparison  of  the  buoy  movement  versus  the  ADCIRC 
currents  (figure  4.12a),  shows  good  agreement  between  the  model  and  the  observations  once 
the  buoy  is  moving  along  its  southeast  track.  When  considering  the  total  currents  (movement 
of  the  buoy  plus  the  ADCP  measured  currents),  the  northwest  currents  leading  up  to  the  peak 
surge  is  well  represented  by  the  model.  Immediately  after  landfall,  the  northward  currents 
from  the  model  track  the  observations  but  the  eastward  component  never  reaches  the 
maximum  observed  currents  and  then  begins  to  diminish  while  the  observed  currents 
continue  to  increase.  Two  potential  causes  can  be  attributed  to  this  discrepancy.  First,  the 
model  is  forced  with  the  NOAA  reanalysis  of  the  core  circulation  of  the  storm  which  does 
not  include  the  far  field  winds  that  continue  offshore  beyond  the  time  of  storm  landfall. 
Hence  the  model  agrees  reasonable  well  with  the  observations  during  the  landfalling  period 
of  the  storm  but  currents  fall  off  rapidly  after  that  point  in  time.  Secondly,  the  model  retained 
very  high  waters  inland  (recall  figure  4.7a)  due  to  shortcomings  in  the  wetting  and  drying 
algorithm  as  exposed  in  section  3.  Thus  the  model  does  not  capture  retreating  waters  as 
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Figure  4. 1 1 :  (a)  Eastward  wind  stress  observed  (blue)  and  applied  to  the  model  (red);  (b) 
Northward  wind  stress  observed  (blue)  and  applied  to  the  model  (red). 


Buoy  Velocity 


Absolute  Currents  (Buoy  +  ADCP) 


Figure  4.12:  Northward  (red)  and  eastward  (green)  movement  of  Buoy  USM3m001  and  the 
model  computed  northward  and  eastward  components  of  the  currents  (magenta  and  cyan, 
respectively);  (b)  the  total  measured  (buoy  movement  plus  ADCP  measured  currents)  and  the 
modeled  currents  for  the  northward  and  eastward  components  (same  color  scheme). 
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evident  by  the  buoy  recorded  currents.  Figure  4.13  provides  maps  of  the  modeled  current 
magnitudes  and  directions  at  two  times  during  the  observed  storm  event  (indicted  by  black 
vertical  lines  on  figure  4.12b. 


Velocity  Magnitude  (m/s)  on  29-August-2005  at  8:00:00  AM 


-90  -89.8  -89.6  -89.4  -89.2  -89  -88.8  -88.6  -88.4  -88.2  -88 

Longitude  (deg)  (a) 


Velocity  Magnitude  (m/s)  on  29- August-2005  at  11:59:59  AM 
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Figure  4.13:  Modeled  current  magnitude  (color)  in  m/s  and  direction  (white  arrows)  at  (a) 

8:00  CDT  29  August  2005  and  (b)  12:00  CDT  29  August  2005.  These  times  are  indicated  by 
black  lines  in  figure  4.13b. 

e.  Summary 

The  ability  to  rapidly  apply  the  ADCIRC  surge  model  to  any  location  globally  is  the 
operational  goal  of  the  surge  forecast  prediction  system.  The  mesh  generation  tool, 
MeshGUI,  renders  that  goal  of  re-locateability  possible.  Experiences  gained  during  the 
hindcast  of  Katrina  have  lead  to  upgrades  in  the  MeshGUI  tool.  For  example,  a  series  of 
mesh  quality  adjustments  are  now  automatically  applied  to  a  created  mesh  to  eliminate 
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poorly  constructed  triangular  elements  that  may  cause  model  instabilities.  Furthermore  an 
estimate  of  model  computational  time  is  provided  based  on  the  size  of  a  created  mesh.  The 
user  can  decide  if  iterations  on  the  mesh  design  are  needed  knowing  current  operational 
constraints.  Additionally,  experience  in  creating  the  “Katrina”  mesh  indicates  that  a  multi¬ 
stage  mesh  generation  approach  may  be  advantageous  to  balance  resolution  requirements  and 
computational  limitations  in  different  regions  of  the  mesh.  For  this  approach  the  ability  to 
“stitch”  together  different  meshed  regions  was  developed.  Our  own  as  well  as  others’ 
experience  modeling  the  inundation  from  Katrina  indicates  that  fine  scale  information  on 
overland  elevation,  vegetation  type,  and  frictional  characteristics  are  all  very  important  for 
accurate  represention  of  a  wetting  front.  Methods  are  now  being  developed  to  automatically 
extract  such  information  from  remotely  sensed  imagery  and  utilize  it  in  the  mesh  generation 
process. 

Apart  from  the  mesh  generation  process,  our  analyses  of  Hurricane  Katrina  surge  and 
inundation  hindcasts  highlight  improvements  to  the  inundation  methodology  that  could  result 
in  even  more  accurate,  robust  forecasts.  For  example,  the  movement  of  water  inland  would 
be  better  represented  as  a  response  to  not  only  water  elevation  and  frictional  effects  but  also 
wind  forcing.  The  conservation  of  water  in  overland  regions  that  are  wet,  dried  and  rewet  is 
another  important  aspect.  In  addition  tracking  and  resolving  the  wet-dry  interface  could 
further  enhance  fidelity  of  the  inundation  forecast.  Resolution  of  the  mesh  and  representation 
of  the  bathymetry  and  topography  as  well  as  aresonable  representation  of  the  wind  event  are 
the  most  essential  elements  for  accurate  storm  surge  and  inundation  prediction. 


5.  Validation  Test  Results  -  Delaware  Bay  Tides  and  Currents 

NOAA’s  National  Ocean  Service  (NOS)  is  in  the  processes  of  evaluating 
oceanographic  nowcast  and  forecast  modeling  systems  to  support  navigational  and 
environmental  applications  in  U.S.  coastal  waters.  In  that  context  NOAA  NOS  has 
established  the  Delaware  Bay  benchmark  for  evaluating  various  model’s  performance  using 
in-house  developed  skill  assessment  software.  We  adopt  the  same  Delaware  Bay  benchmark 
for  validation  of  the  3D  barotropic  ADCIRC  model  in  the  spirit  of  using  community 
established  benchmarks  and  to  take  advantage  of  the  rich  data  set  associated  with  this 
benchmark.  The  objectives  of  the  Delaware  Bay  benchmark  as  applied  here  are  1)  to  validate 
and  verify  ADCIRC  elevation  and  2D  and  3D  current  predictions  against  observational  data 
and  2)  to  evaluate  the  model’s  response  and  associated  uncertainty  in  an  estuary  setting  under 
various  forcing  scenarios. 

Delaware  Bay  is  a  major  estuary  of  the  U.S.  east  coast,  surrounded  by  the  states  of 
Delaware,  Pennsylvania  and  New  Jersey.  The  length  of  the  bay  is  about  75  km,  has  a 
maximum  width  of  about  45  km  with  a  mean  depth  of  10  m  (figure  5.1a).  The  study  area 
spans  the  distance  from  the  Atlantic  Ocean  to  Trenton,  New  Jersey.  Water  motions  in  the  bay 
are  dominated  by  tidal  currents  with  M2  the  dominant  tidal  constituent.  Delaware  Bay  adds  a 
tidally-driven,  shallow  estuary  with  river  contribution  to  the  suite  of  validation  benchmarks. 

a.  Model  Configuration 

The  unstructured  finite  element  grid  shown  in  figure  5.1b  consists  of  15726  nodes  and 
28831  elements  with  resolution  ranging  from  50  meters  at  the  upstream  tributary  area  to  1000 


42 


meters  at  the  offshore  open  boundary  region.  The  domain  includes  the  Delaware  River,  the 
Chesapeake  and  Delaware  (C&D)  Canal  and  extends  out  from  the  estuary  into  the  Atlantic 
Ocean  to  depths  of  100m.  Forcing  is  applied  in  the  form  of  13  tidal  harmonics  (M2,  S2,  N2, 
Ki,  Oi,  K2,  Qi,  M4,  Mg,  MN4,  Pi,  and  2SM2)  along  the  model’s  open  water  boundary  and  at 
the  C&D  Canal.  Values  for  the  harmonic  constants  are  taken  from  the  2001  ADCIRC  East 
Coast  database  (EC2001)  (Feyen  and  Yang,  2007).  River  discharge  specified  through 
ADCIRC ’s  fort. 20  file  forced  the  head  of  the  Delaware  River.  For  most  of  the  validation 
testing  presented,  the  affects  of  meteorological  forcing  are  ignored.  Two  short  term 
simulations  demonstrate  the  influence  winds  for  the  two  dominant  wind  patterns,  from  the 
northwest  and  from  the  southeast. 

The  ADCIRC  model  was  run  for  190  day  period  in  1984,  a  time  period  during  which 
17  NOS  water  level  gauges  (figure  5.2a)  and  44  current  stations  (figure  5.2b)  recorded 
observations.  A  range  of  experiments  were  configured  to  examine  different  aspects  of  the 
model  forcing  (river  discharge  and  wind)  and  bottom  friction  coefficient  specification.  For 
this  later  set  of  experiments  used  a  spatially  varying  bottom  friction  coefficient  as  suggested 
by  Walters  (1997)  such  that  a  value  of  0.0025  was  specified  in  the  bay  and  the  friction 
increased  to  0.04  upstream  in  the  Delaware  River.  Table  1  presents  a  complete  list  of  the 
numerical  experiments.  Experiments  1  to  6  focused  on  two-dimensional  depth-integrated 
solutions  from  ADCIRC  while  experiments  7  to  12  were  analyzing  three-dimensional 
calculations.  Standard  statistical  measurements  were  computed  to  evaluate  model  skill.  The 
metrics  include  Mean  Bias  Error  (MBE),  Standard  Deviation  (SD),  Root  Mean  Square  Error 
(RSME)  and  Index  of  Agreement  (AI)  as  defined  by  Blain  (1997). 


Table  5. 1  Description  of  the  2D  and  3D  validation  experiments  for  Delaware  Bay. 


Experiment 

No. 

Dimensionality 

Description 

Label 

1 

2D 

Baseline  (constant  river  inflow) 

2DCR 

2 

2D 

No  river  inflow 

2DNR 

3 

2D 

Time-varying  real  river  inflow 

2DRR 

4 

2D 

Spatially  varying  bottom  friction 

2DBF 

5 

2D 

Dominant  NW  wind  event  at  12 m/s 

2DNW 

6 

2D 

Dominant  SE  wind  event  at  12 m/s 

2DSE 

7 

3D 

Baseline  (constant  river  inflow) 

3DCR 

8 

3D 

No  river  inflow 

3DNR 

9 

3D 

Time -varying  real  river  inflow 

3DRR 

10 

3D 

Spatially  varying  bottom  friction 

3DBF 

11 

3D 

Dominant  NW  wind  event  at  12  m/s 

3DNW 

12 

3D 

Dominant  SE  wind  event  at  12  m/s 

3DSE 
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Bathymetry  (+)  /  Topography  (-) 


Figure  5.1:  (a)  Bathymetry  contours  in  m  for  Delaware  Bay  and  (b)  the  unstructured  finite 
element  mesh  created  for  the  same  region  as  part  of  the  NOAA’s  NOS  Delaware  Bay 
benchmark  test  case. 


ADCIRC  Delaware  Bay  Validation,  Tidal  Elevation  Stations  (17)  ADCIRC  Delaware  Bay  Validation,  Tidal  Current  Stations  (44) 


Figure  5.2:  Locations  of  (a)  17  tide  gauges  and  (b)  44  current  meters  associated  with  the 
NOAA  Delaware  Bay  benchmark  test  case. 
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b.  2D  Harmonic  Tidal  Elevations 


Model-data  comparisons  for  the  baseline  NOAA  case  (2D,  constant  river  inflow, 
2DCR)  resulted  in  a  small  0.018m  error  for  the  mean  bias  error  (MBE)  and  a  very  high 
agreement  index  (AI)  equal  to  0.94.  The  influences  of  river  discharge  and  bottom  friction  on 
predictive  skill  are  then  investigated  for  three  different  river  inflow  conditions  (No  River, 
Constant  River,  and  Time- varying  Real  River)  and  variable  bottom  friction  coefficients. 
Figure  5.3  and  figure  5.4  compare  four  error  measures,  mean  bias  error  (MBE),  standard 
deviation  (SD),  mean  absolute  error  (MAE),  and  root  mean  square  error  (RMSE),  for  each  of 
eight  elevation  harmonic  constituents  in  amplitude  and  phase,  respectively.  Overall,  the 
dominate  tidal  constituent,  M2,  has  errors  of  about  6cm,  an  increase  over  the  constant  river 
inflow  condition  Errors  of  about  2 cm  are  computed  for  the  S2,  N2,  Ki  and  M4  constituents 
and  approximately  1cm  error  for  Oi,  M6  and  MN4  constituents  relative  to  the  2DCR  and 
2DNR  cases.  Changing  the  bottom  friction  coefficient  from  0.0025  to  a  spatially  varying 
bottom  stress  produces  the  largest  errors.  For  the  M2  constituent,  MBE  =-  0.20m,  SD  =  0.23 
m,  MAE  =  0.22  m,  RMSE  =  0.30  m  and  AI=  0.35.  All  four  cases  (2DCR,  2DNR,  2DRR, 
2DBF)  have  similar  magnitudes  for  the  phase  error  as  seen  in  figure  5.4.  There  is  a  10  degree 
variation  for  M2,  S2,  N2,  and  Oi  and  much  more  substantial  phase  errors  in  the  range  of  40 
degrees  for  M4,  M6,  and  MN4,  the  shallow  water  quarter  diurnal  constituents. 

The  amplitude  Agreement  Index  ranges  from  0.8  to  0.95  for  all  constituents  in  each  of 
the  test  cases  having  variants  of  the  river  forcing  with  the  exception  of  the  Ki  constituent. 
Varying  the  bottom  stress  resulted  in  notably  lower  AI  scores.  For  phase,  all  cases  tested 
have  excellent  phase  agreement  with  observations  according  to  the  AI  value  which  ranges 
from  0.95  to  0.99  for  dominant  as  well  as  secondary  constituents  (figure  5.5). 

Overall,  whether  no  river  discharge  or  a  constant  river  discharge  are  specified  all 
harmonic  tidal  elevation  errors  are  quite  similar.  Adding  a  time -varying  daily  discharge 
which  is  more  representative  of  reality  did  lead  to  slight  improvements  in  the  model-data 
comparisons.  Spatially  varying  the  bottom  friction,  as  suggested  by  Walters  (1997)  to 
improve  diurnal  constituent  representation,  leads  to  large  error  with  a  negative  bias  indicating 
the  model  computations  are  very  sensitive  to  the  bottom  friction  coefficient.  Increasing 
bottom  friction  heavily  damped  the  elevation  signal  and  degraded  the  comparison  to 
observed  water  levels. 

To  investigate  any  spatial  trends  in  the  water  level  error,  individual  station  errors  are 
computed  considering  only  the  dominant  tidal  constituent,  M2.  Table  5.2  summarizes  the 
individual  station  errors  for  each  experiment.  These  same  errors  are  graphically  depicted  in 
figure  5.6.The  17  stations  are  grouped  into  three  categories  depending  on  their  location  to  see 
if  any  trend  or  bias  can  be  found.  Four  stations  are  in  the  estuarine  category  (8539993, 
8539487,  8539094,  8545530)  indicating  primary  influence  of  the  Delaware  River,  eight 
stations  are  in  the  bay  category  (8551762,  8551910,  8573927,  8537614,  8554399,  8555388, 
8536581,  85555889),  and  five  stations  (8534720,  8536110,  8557380,  8558690,  8570280)  are 
identified  as  offshore  near  the  open  ocean  boundary.  For  the  constant  river  flow  case,  the 
model  under-predicts  water  level  at  estuarine  and  bay  stations  while  over-estimating  water 
height  at  open  ocean  stations.  But  for  the  2DNR  and  2DRR  experiments,  the  model  over¬ 
predicts  at  estuarine  and  open  ocean  stations  and  still  slightly  under-predicts  water  level  at 
the  bay  stations.  For  the  2DBF  case,  the  model  under-estimated  water  levels  at  all  but  two 
stations. 

Among  all  17  stations,  the  largest  magnitude  errors  are  located  at  Trenton  (8539993), 
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Artificial  Island  (8537614)  and  Indian  River  Inlet  (8558690).  Errors  at  these  three  stations 
are  at  least  one  order  of  magnitude  larger  than  errors  at  the  remaining  stations.  The  probable 
cause  for  the  large  error  at  Indian  River  Inlet  (found  outside  the  Bay  along  the  western 
shoreline)  is  the  complex  inlet  geometry  not  currently  included  in  the  model  grid. 
Additionally,  discharge  from  the  Indian  River  is  not  incorporated  as  forcing.  At  the  Trenton 
station  (located  at  the  northeastern  most  point  of  the  Delaware  River)  error  is  due  to  the 
variability  in  the  Delaware  river  flux.  When  a  time-varying  river  flux  is  applied  as  forcing, 
the  error  reduces  from  10cm  down  to  1cm.  Similar  errors  are  computed  for  the  constant  and 
no  river  flux  forcing  cases.  Errors  at  Artificial  Island  located  just  south  of  the  Chesapeake  & 
Delaware  Canal  are  due  to  combined  errors  in  the  forcing  specified  for  both  the  Delaware 
River  and  the  C&D  canal.  A  spatially  varying  bottom  friction  still  creates  the  largest  errors  at 
nearly  every  station,  especially  in  the  estuarine  and  bay  regions  where  bottom  friction 
coefficients  are  higher. 

Table  5.2  2D  tidal  elevation  errors  for  the  M2  tide  at  17  water  level  stations  for  all 


experiments. 


Location 

StationID 

2DCR 

2DNR 

2DRR 

2DBF 

Trenton* 

8539993 

-0.1065 

0.0687 

0.0150 

-0.8595 

Fieldsboro 

8539487 

-0.0066 

0.1484 

0.1092 

-0.7280 

Burlington 

8539094 

-0.0532 

0.0694 

0.0563 

-0.6832 

Philadelphia 

8545530 

-0.0483 

0.0407 

0.0441 

-0.2364 

Delaware  City 

8551762 

-0.0629 

-0.0372 

-0.038 

-0.1796 

Reedy  Pt. 

8551910 

-0.0784 

-0.0599 

-0.0623 

-0.1880 

Chesapeake  City 

8573927 

-0.0090 

-0.0132 

-0.0142 

-0.0039 

Artificial  Island 

8537614 

-0.1654 

-0.1354 

-0.1377 

-0.2805 

Mahon  Riv 

8554399 

0.0007 

0.0360 

0.0384 

-0.1500 

Murderkill  Riv 

8555388 

0.0359 

0.0702 

0.0695 

-0.1235 

Bidwell  Creek 

8536581 

-0.0013 

0.0349 

0.0346 

-0.1671 

Brandywine 

8555889 

-0.0003 

0.0276 

0.0274 

-0.1326 

Atlantic  City 

8534720 

0.0036 

0.0036 

0.0036 

0.0036 

Cape  May  Canal 

8536110 

0.0069 

0.0332 

0.0331 

-0.1120 

Lewes 

8557380 

0.0153 

0.0328 

0.0327 

-0.066 

Indian  Riv  Inlet 

8558690 

0.1825 

0.1852 

0.1852 

0.1674 

Ocean  City 

8570280 

0.0105 

0.0106 

0.0106 

0.0105 
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Figure  5.3:  Four  2D  amplitude  error  measures,  mean  bias  error  (MBE),  standard  deviation 
(SD),  mean  algebraic  error  (MAE),  and  root  mean  square  error  (RMSE),  for  each  of  8  water 
level  harmonic  tidal  constituents  for  the  (a)  2DCR,  (b)  2DNR,  (c)  2DRR,  and  (d)  2DBF 
experiments. 
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Figure  5.4:  Four  2D  phase  error  measures,  mean  bias  error  (MBE),  standard  deviation  (SD), 
mean  algebraic  error  (MAE),  and  root  mean  square  error  (RMSE),  for  each  of  8  water  level 
harmonic  tidal  constituents  for  the  (a)  2DCR,  (b)  2DNR,  (c)  2DRR,  and  (d)  2DBF 
experiments. 
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dgure  5.5:  The  2D  Index  of  Agreement  (AI)  elevation  (a)  amplitude  and  (b)  phase  for  each 
of  8  water  level  harmonic  tidal  constituents  for  the  2DCR,  2DNR,  2DRR,  and  2DBF 
experiments. 
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Error  for  2D  no  river  case 


Error  for  2D  const,  river  case 


Error  for  2D  real  river  case  Error  for  2D  friction  case 


Figure  5.6:  The  M2  elevation  MBE  for  amplitude  from  each  of  the  2D  experiments  (a) 
2DCR,  (b)  2DNR,  (c)  2DRR,  and  (d)  2DBF. 
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c.  3D  Harmonic  Tidal  Elevations 


A  similar  analysis  as  outlined  in  section  5b  was  conducted  using  results  from  the  3D 
ADCIRC  computations.  Amplitude  and  phase  errors  for  8  constituents  over  4  3D  model 
experiments  are  shown  in  figures  5.7  and  5.8.  For  the  3D  experiments,  errors  increase  for  the 
M2  (10 cm)  and  M4  (5 cm)  constituents  when  compared  with  the  2D  results.  The  baseline 
(2DCR)  experiment,  had  errors  for  the  M2  tidal  constituent  as  follows:  a  MBE  of  -0.0028m, 
SD  of  0.104m,  MBE  of  0.0812m,  an  RMSE  of  0.102m  and  an  AI  value  of  0.849.  The  lowest 
error  values  are  recorded  for  the  3DRR  experiment  for  the  M2  tidal  constituent  when 
compared  to  the  other  experiments  and  follows  the  same  trend  as  the  2D  error  analysis.  In 
contrast  notably  larger  errors  are  associated  with  the  S2  amplitude  for  the  3DRR  test  case. 
Other  constituents  (N2,  Ki,  M4,  Oi,  Mg  and  MN4)  exhibit  similar  statistical  scores  for  all  test 
cases  though  errors  are  increased  for  the  higher  harmonics  (M4,  M6  and  MN4).  Phase  errors 
in  all  experiments  are  under  10  degrees  with  the  exception  of  the  higher  harmonics  (M4,  M6 
and  MN4).  No  significant  differences  in  the  phase  error  are  evident  among  the  test  cases.  The 
Agreement  Index  for  amplitude  and  phase  for  the  3D  experiments  are  provided  in  figure  5.9. 
Values  of  the  AI  are  quite  similar  to  those  computed  from  the  2D  model  data  except  for  the 
quarter-diurnal  constituents,  M4  and  MN4.  All  constituents  in  3D  runs  also  show  excellent 
phase  agreement  with  AI  ranges  from  0.9  to  0.99. 

Station  by  station  error  values  for  17  water  level  stations  are  listed  in  Table  5.3  for  the 
M2  tide  for  all  experiments.  For  3D  experiments  the  largest  errors  occur  at  the  estuarine 
stations  (Trenton,  Fieldsboro  and  Burlington)  for  the  3DCR,  3DNR,  and  3DBF  test  cases. 
Error  magnitudes  for  the  3DCR  and  3DNR  are  similar  at  most  stations.  Significant 
improvement  is  evident  at  the  4  estuarine  stations  for  the  3DRR  but  no  improvement  is  seen 
at  either  the  bay  or  open  ocean  stations.  When  comparing  differences  between  the  2D  and  3D 
constant  river  flow  case,  the  2DCR  actually  has  smaller  errors  at  the  estuarine  stations  than 
the  3D  case.  From  this  set  of  experiments,  the  3D  model  performs  best  for  the  river  forcing 
that  represents  the  more  realistic  time-varying  river  discharge. 

Spatially,  the  3DCR,  3DNR,  and  the  3DBF  experiments  show  under-prediction  at  the 
estuarine  stations  and  over-prediction  at  the  open  ocean  stations.  The  exception  is  the  3DRR 
which  has  a  positive  bias  for  all  estuarine  stations. 
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Table  5.3  3D  tidal  elevation  errors  for  the  M2  tide  at  17  water  level  stations  for  all 


experiments. 


Location 

StationID 

3DCR 

3DNR 

3DRR 

3DBF 

Trenton* 

8539993 

-0.2502 

-0.2546 

0.0180 

-0.2502 

Fieldsboro 

8539487 

-0.1695 

-0.1746 

0.0934 

-0.1695 

Burlington 

8539094 

-0.1931 

-0.2124 

0.0394 

-0.1931 

Philadelphia 

8545530 

-0.090 

-0.1409 

0.0621 

-0.090 

Delaware  City 

8551762 

-0.0013 

0.0075 

0.0271 

-0.0013 

Reedy  Pt. 

8551910 

-0.0158 

-0.0082 

-0.0066 

-0.0158 

Chesapeake  City 

8573927 

-0.010 

-0.011 

-0.0090 

-0.0102 

Artificial  Island 

8537614 

-0.1067 

0.0986 

-0.0916 

-0.1067 

Mahon  Riv 

8554399 

0.0646 

0.0673 

0.0827 

0.0646 

Murderkill  Riv 

8555388 

0.1109 

0.1124 

0.1278 

0.1109 

Bidwell  Creek 

8536581 

0.0723 

0.0731 

0.0875 

0.0723 

Brandywine 

8555889 

0.0718 

0.0724 

0.0815 

0.0718 

Atlantic  City 

8534720 

0.0036 

0.0036 

0.0036 

0.0036 

Cape  May  Canal 

8536110 

0.0666 

0.0670 

0.0812 

0.0666 

Lewes 

8557380 

0.0512 

0.0513 

0.0580 

0.0512 

Indian  Riv  Inlet 

8558690 

0.1795 

0.1794 

0.1802 

0.1795 

Ocean  City 

8570280 

0.0104 

0.0104 

0.0104 

0.0104 
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Figure  5.7:  Four  3D  amplitude  error  measures,  mean  bias  error  (MBE),  standard  deviation 
(SD),  mean  algebraic  error  (MAE),  and  root  mean  square  error  (RMSE),  for  each  of  8  water 
level  harmonic  tidal  constituents  for  the  (a)  3DCR,  (b)  3DNR,  (c)  3DRR,  and  (d)  3DBF 
experiments. 
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Figure  5.8:  Four  3D  phase  error  measures,  mean  bias  error  (MBE),  standard  deviation  (SD), 
mean  algebraic  error  (MAE),  and  root  mean  square  error  (RMSE),  for  each  of  8  water  level 
harmonic  tidal  constituents  for  the  (a)  3DCR,  (b)  3DNR,  (c)  3DRR,  and  (d)  3DBF 
experiments. 
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Figure  5.9:  The  3D  Index  of  Agreement  (AI)  for  elevation  (a)  amplitude  and  (b)  phase  for 
each  of  8  water  level  harmonic  tidal  constituents  for  the  3DCR,  3DNR,  3DRR,  and  3DBF 
experiments. 
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Error  for  3D  no  river  case 


Error  for  3D  const,  river  case 


Error  for  3D  real  river  case  Error  for  3D  friction  case 


Figure  5.10:  The  M2  elevation  MBE  for  amplitude  from  each  of  the  3D  experiments  (a) 
3DCR,  (b)  3DNR,  (c)  3DRR,  and  (d)  3DBF. 
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d.  2D  Harmonic  Tidal  Currents 


Statistical  analyses  of  the  harmonic  constituent  amplitudes  and  phases  of  2D  currents 
are  presented  for  the  four  experiments,  2DCF,  2DNR,  2DRR,  and  2DBF.  The  constant  river 
flow  case  (2DCR)  is  used  as  a  benchmark  and  is  compared  against  the  other  three  test  cases. 
The  average  errors  are  20cm/s  for  M2  and  5 cm/s  for  the  remaining  constituents.  Specification 
of  a  spatially  variable  bottom  friction  once  again  yields  larger  MBE  and  RMSE  values  as  it 
over-damps  the  velocity  in  the  water  column.  All  current  tidal  constituents  for  each  of  the 
four  test  cases  had  phase  errors  ranging  from  70  to  120  degrees  with  no  significant 
differences  regardless  of  the  river  forcing  or  the  spatially  variable  bottom  friction 
specification.  Some  of  the  phase  error  may  be  attributed  to  inconsistencies  in  the  axis 
orientation  used  for  comparing  observations  and  modeled  results.  (The  dominant  M2  tide 
enters  the  estuary  at  110  degree  from  true  north,  with  a  clockwise  rotation.  Observations  are 
processed  using  the  true  north/east  directions  while  the  skill  assessment  software  may 
decompose  the  currents  to  major/minor  axis.  A  110  degree  phase  shift  in  the  modeled 
major/minor  axis  results  in  phases  errors  of  20-30  degree  which  is  more  in  keeping  with  the 
known  model  skill. 

The  temporally  varying  real  river  simulation  provided  the  best  overall  statistical  skill. 
For  2DRR,  the  M2  amplitude  had  the  smallest  MBE;  it  also  had  the  smallest  SD,  MAE, 
RSME  and  highest  AI  among  all  four  cases  for  most  of  the  tidal  constituents.  2DCR  and 
2DNR  yielded  comparable  statistical  numbers  while  the  varying  bottom  case  (2DBF) 
generated  larger  MBE,  RSME  and  lowest  AI  for  the  M2  amplitude.  Figure  5.1 1  shows  the  AI 
score  ranging  from  0.6  to  0.8  for  current  amplitude  and  phase. 

The  spatial  distribution  of  M2  current  errors  are  presented  in  figure  5.12.  Unlike  the 
water  level  errors  discussed  in  the  previous  sections,  large  errors  for  the  currents  are 
concentrated  in  the  upper  bay  where  maximum  velocities  usually  occur  in  conjunction  with 
complex  stratification  and  mixing  processes.  For  2DCR,  open  ocean  stations  have  the 
smallest  error,  between  3 cm  and  5 cm,  estuarine  stations  have  errors  between  5  and  15cm 
while  errors  at  bay  stations  are  as  high  as  60cm.  By  incorporating  the  realistic  time -varying 
river  flow,  velocity  errors  are  reduced  in  both  the  estuarine  and  bay  regions.  No  improvement 
is  observed  for  the  2DBF  case,  despite  the  increased  bottom  friction  in  the  estuarine  and  bay 
regions. 
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Figure  5.11:  The  3D  Index  of  Agreement  (AI)  for  current  (a)  amplitude  and  (b)  phase  for 
each  of  8  water  level  harmonic  tidal  constituents  for  the  3DCR,  3DNR,  3DRR,  and  3DBF 
experiments. 
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Figure  5.12:  The  M2  velocity  MBE  for  each  of  the  2D  experiments  (a)  2DCR,  (b)  2DNR,  (c) 
2DRR,  and  (d)  2DBF. 
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e.  3D  Harmonic  Tidal  Currents 


Average  errors  for  the  model  representation  of  3D  currents  in  the  constant  river  case 
are  20 cm/s  for  M2  and  5 cm/s  for  the  rest  of  constituents,  the  same  magnitude  of  error  as  for 
the  2D  simulation.  Variable  bottom  friction,  3DBF,  results  in  larger  MBE  and  RMSE  values. 
All  constituents  in  each  of  the  four  cases  have  phase  errors  ranging  from  100  to  150  degrees 
with  no  significant  differences  observed  between  the  various  experiments.  Once  again,  the 
realistic  river  inflow  (3DRR)  demonstrated  the  best  model  skill.  The  AI  values  for  current 
amplitudes  from  the  3DRR  experiment  are  higher  than  the  remaining  three  test  cases  in  9  of 
13  constituents,  including  the  major  M2,  S2,  N2  and  Ki  constituents  (see  figure  5.13).  The  AI 
for  current  phase  for  the  3DRR  experiment  is  also  higher  for  8  of  13  constituents  when 
compared  with  other  three  test  cases.  3DRR  also  has  the  smallest  SD,  MAE,  RSME  and 
highest  AI  among  all  four  cases  for  most  of  the  tidal  constituents.  The  3DCR  and  3DNR 
cases  yield  comparable  statistical  errors  while  the  varying  bottom  case  (3DBF)  still  results  in 
larger  MBE  and  RSME  and  the  lowest  AI  values.  AI  values  for  current  amplitudes  range 
from  0.6  to  0.8  while  AI  for  current  phases  have  a  wider  range  from  0.4  to  0.8  in  figure  5.13. 
The  error  analyses  have  shown  the  importance  of  the  inclusion  of  realistic  river  discharge  in 
both  the  2D  and  3D  model  computations,  even  model  computations  neglecting  density 
contributions. 

Spatially  the  errors  for  the  3DCR  test  case  showed  a  reduction  in  the  model-data  error 
in  the  upper  bay  region,  but  not  in  the  estuarine  section  when  compared  with  errors  computed 
using  the  2D  solution  for  the  same  forcing  (see  figure  5.14).  The  3DRR  case  had  the  smallest 
errors  in  both  the  estuarine  and  bay  regions  as  compared  with  the  2DCR,  2DRR  and  the 
3DCR  experiments. 
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Figure  5.13:  The  3D  Index  of  Agreement  (AI)  for  current  (a)  amplitude  and  (b)  phase  for 
each  of  8  water  level  harmonic  tidal  constituents  for  the  3DCR,  3DNR,  3DRR,  and  3DBF 
experiments. 
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M2  Velocity  Error  for  3D  CR  M2  Velocity  Error  for  3D  NR 


Figure  5.14:  The  M2  velocity  MBE  for  each  of  the  3D  experiments  (a)  3DCR,  (b)  3DNR,  (c) 
3DRR,  and  (d)  3DBF. 
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f.  Water  Level  Time  Series  Comparisons 


Of  all  the  17  NOAA  water  level  stations,  only  seven  are  primary  stations  with  long¬ 
term  continuous  records.  There  are  two  located  on  the  coast,  two  within  Delaware  Bay  and 
three  near  the  Delaware  River.  These  seven  NOAA  stations  with  complete  water  level  data 
records  are  used  for  model-data  time  series  comparisons.  Twenty  five  days  of  simulated 
water  levels  for  each  of  the  four  test  cases  (2DCR,  2DNR,  2DRR,  and  2DBF)  are  compared 
with  the  observation  data.  Mean  model  and  observed  water  levels  are  computed  and  mean 
algebraic  errors  are  reported  in  Table  5.4.  In  general,  ADCIRC  generated  water  levels  are  in 
good  agreement  with  measured  data,  for  both  amplitude  and  phase.  The  largest  errors  occur 
at  two  river  inlet  locations.  The  form  of  the  river  forcing  (2DCR,  2DNR,  or  2DRR)  does  not 
make  a  significant  difference  in  the  water  level  time  series,  although  minor  fluctuations  are 
evident  nearest  the  river  discharge  location.  The  2D  and  3D  modeled  water  levels  simulations 
produced  similar  results.  The  most  notable  influence  came  from  the  applied  meteorological 
forcing.  Strong  wind  events  (12 m/s)  showed  a  -25  to  +25  cm  variation  when  compared  to 
model  computations  without  surface  wind.  The  magnitude  of  the  deviation  depends  on  the 
wind  speed  and  direction. 


Table  5.4.  Mean  algebraic  error 

(MAE)  for  water  level  at  7  NOAA  stations  for  all  test  cases. 

Station 

2DCR 

2DRR 

2DNR 

2DBF 

(MAE) 

(MAE) 

(MAE) 

(MAE) 

8534720 

0.0835 

0.0834 

0.0834 

0.0834 

8558690 

0.0379 

0.0386 

0.0385 

0.0390 

8536110 

0.0245 

0.0239 

0.0237 

0.0358 

8557380 

0.0273 

0.0296 

0.0294 

0.0311 

8551910 

0.1074 

0.1423 

0.1369 

0.1313 

8545240 

0.1607 

0.2131 

0.1854 

0.2146 

8539993 

0.2346 

0.1831 

0.1129 

0.9859 

The  model-data  time  series  are  plotted  for  the  2DCR  baseline  at  all  7  NOAA  stations 
in  figure  5.15  for  coastal  and  bay  stations  and  figure  5.16  for  river  stations.  At  coastal  station 
8534720  (figure  5.15a),  the  model  represents  well  the  observed  water  levels  after  a  ramping 
period  of  5  days.  In  this  location,  model  predictions  and  observational  data  are  in  good 
agreement  with  an  MAE  value  of  just  above  8 cm.  No  notable  phase  lag  or  lead  is  evident  as 
the  modeled  water  levels  track  the  observed  peaks  and  troughs.  Errors  appear  larger  at  the 
other  coastal  station  (8558690)  shown  in  figure  5.15b  thought  he  MAE  has  a  smaller  value  of 
0.0379m.  The  model  overestimates  water  levels  much  of  the  time  series  with  observed  water 
levels  having  a  smaller  tidal  range.  Model  computations  at  two  bay  stations  (8536110  and 
8557380),  shown  in  figure  5.15c  and  5.15d,  have  the  lowest  MAE  just  above  2 cm.  At  these 
stations  the  model  underestimates  high  tide  during  the  first  eight  days,  then  overestimates  the 
low  tide  elevation  for  the  next  few  days.  For  the  final  5  days  the  model  continues  to 
underestimate  the  low  tide  but  grossly  overestimates  the  high  tide  during  the  same  period. 
Tides  at  the  three  river  stations  had  stronger  tidal  signals  when  compared  to  water  level  time 
series  at  the  other  four  stations;  MAE  values  range  from  10cm  to  23 cm.  Station  8551910 
showed  reasonable  agreement  between  model  and  measured  data  (figure  5.16a)  and  follows 
as  error  pattern  similar  to  bay  station  8536110  (figure  5.15c).  In  contrast,  at  stations  8539993 
and  8545530,  the  model  computed  water  levels  overestimate  high  tide  and  underestimate  low 
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tide  (figures  5.16b  and  5.16c). 


Time  Series  comparison  at  N  0$  station  8534720 


Time  Series  comparison  at  NOS  station 85361 10 
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Figure  5.15:  Water  level  times  series  of  model  (black)  versus  observations  (red)  for  the 
2DCR  experiment  at  4  NOAA  stations,  on  the  coast  (a)  8534720,  (b)  8558690,  and  in  the  bay 
(c)  8536110,  and  (d)  8557380. 
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Time  Series  comparison  at  NOS  station  8551910 


Figure  5.16:  Water  level  times  series  of  model  (black)  versus  observations  (red)  for  the 
2DCR  experiment  at  3  NOAA  stations  near  the  river  (a)  8551910,  (b)  8545240,  and  (c) 
8539993. 

Comparisons  between  the  40-day  time  series  for  all  four  experiments  (2DCR,  2DNR, 
2DRR,  and  2DBF)  show  negligible  differences  for  the  various  river  discharge  options.  The 
water  levels  are  almost  identical  at  open  water  stations  where  river  flux  is  negligible.  Even  at 
a  location  near  the  Delaware  River  mouth,  the  difference  between  no  river  discharge  and 
realistic  river  inflow  is  only  few  centimeters  (figure  5.18e).  As  noted  in  prior  discussions, 
modeled  elevations  are  extremely  sensitive  to  the  bottom  friction  coefficient.  The  increased 
bottom  friction  coefficient  in  the  estuarine  region  damps  the  tidal  signal  to  a  range  of  0.5m 
(figure  5.17e),  a  significant  reduction  in  the  2.5 m  range  of  water  level  fluctuation  for  the 
baseline,  2DCR  test  case. 
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Time  Series  comparison  for  4  eases  at  NOS  station  85361 10 


Time  Series  comparison  for  4  eases  at  N(  >S  station  8557380 


Time  Series  comparison  for  4  coses  at  NOS  Elation  5545530 
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Time  Series  comparison  for  4  eases  at  NOS  station  8551910 
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Time  Series  comparison  for  4  eases  at  NOS  station  8539993 
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Figure  5.17  Modeled  water  level  times  series  for  the  2DCR  (black),  2DNR  (red),  2DRR 
(blue),  and  2DBF  (green)  experiments  at  five  NOAA  stations  near  the  river  (a)  85361 10,  (b) 
8557380,  (c)  8545530,  (d)  8551910,  and  (e)  8539993. 
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To  determine  the  influence  of  meteorological  forcing,  strong  wind  events  based  on 
two  predominant  wind  direction  patterns  are  selected  for  simulation.  Southeast  and  northwest 
winds  at  12 m/s  are  applied  to  a  30  day  model  simulation.  The  modeled  elevations  are  not 
compared  to  field  data  since  winds  are  idealized.  Time  series  of  the  computed  elevations  are 
plotted  against  the  no  wind  condition  (case  2DCR)  to  demonstrate  the  effect  of  winds.  As  we 
might  expect,  the  water  levels  generated  without  wind  stress  should  fall  in  magnitude 
between  water  levels  subject  to  winds  in  the  opposite  directions.  Under  the  SW  winds,  a 
strong  water  elevation  setup  occurs  at  two  river  stations  (8539993  and  8545530)  while  water 
levels  at  the  coastal  stations  do  not  exhibit  much  response  to  the  wind  events.  Time  series  for 
the  SE  wind,  the  NW  wind  and  no  wind  are  shown  in  figure  5.18.  The  applied  winds  lead  to  a 
50cm  range  in  the  amplitude  variation  as  compared  to  a  scenario  without  wind.  The 
magnitude  of  this  deviation  depends  on  the  wind  speed  and  direction. 


Wind  impact  on  dcvation  at  NOS  station  8545530 


Wind  impact  on  dcvation  at  NOSstation  8551910 


_i  5 1 _ I _ I _ I _ I _ I _ I _ I _ I _ I _ 

0  2  4  6  8  10  12  14  16  18  20 


Figure  5.18  Modeled  water  levels  at  station  (a)  8545530,  (b)  8551910,  and  (c)  8539993  for 
an  applied  SW  wind  (blue),  NW  wind  (red),  or  no  wind  (black)  for  the  2DCR  test  case. 
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g.  Velocity  Time  Series  Comparisons 


The  evaluation  of  modeled  currents  involves  comparisons  to  current  meter 
observations  at  8  locations.  Three  stations  are  located  at  the  bay  entrance  (001,  003  and  005), 
two  are  inside  the  bay  (023  and  033)  and  three  stations  are  in  the  Delaware  River  (052,  053 
and  054).  A  twenty-day  time  series  of  velocity  magnitudes  are  used  to  compute  average 
velocity  values  at  each  location  as  well  as  a  MAE.  Tabulated  results  are  recorded  in  table  5.5 
for  all  four  test  cases,  2DCR,  2DNR,  2DRR,  and  2DBF.  Time  series  of  the  current 
magnitude  compared  during  periods  where  current  observations  are  available,  e.g.  for 
stations  001,  003,  005,  023,  033,  Julian  Days  95  to  115  are  compared  and  for  stations  052, 
053  and  054,  Julian  Days  67  to  87  are  considered. 

The  average  velocity  for  all  eight  current  meters  ranges  from  0.42  to  0.67  m/s  except 
at  station  054  in  the  upper  Delaware  river  where  currents  are  quite  small  at  0.15mA.  Even 
under  the  constant  river  discharge  (2DCR),  the  modeled  currents  agree  fairly  well  with  the 
observations  having  an  error  of  less  than  5 cm/s  at  5  of  8  locations.  The  largest  discrepancies 
occur  at  station  001,  033  and  053.  For  realistic  river  discharge  (2DRR),  the  model  computed 
currents  improve  compared  to  the  observations  and  errors  are  reduced  at  all  stations,  not  just 
the  upstream  river  stations.  It  should  also  be  noted  that  errors  indicate  that  the  model 
underestimates  current  velocity  at  6  of  the  8  locations.  Not  surprisingly,  modeled  currents 
under  no  river  discharge  as  compared  to  the  realistic  discharge  test  case  yielded  comparable 
results  except  at  the  two  upstream  locations  where  the  2DNR  case  had  larger  errors.  Just  as 
with  water  levels,  the  velocity  signals  are  damped  by  the  increased  bottom  friction  in  the  bay 
and  river  regions  resulting  in  significant  reductions  in  the  velocity  magnitude  and  large  errors 
at  all  stations.  Time  series  of  the  velocity  magnitude  for  the  2DCR  and  2DRR  test  cases  at 
three  stations  are  shown  in  figures  5.19  through  5.21. 


Table  5.5.  Average  velocities  for  modeled  and  observed  currents  and  mean  algebraic  error 
(MAE)  at  8  NOAA  current  meter  sites  for  all  test  cases  (2DCR,  2DNR,  2DRR,  and  2DBF). 


Stn 

Vavg 

(m/s) 

2DCR 

2DCR 

(MAE) 

2DRR 

2DRR 

(MAE) 

2DNR 

2DNR 

(MAE) 

2DBF 

2DBF 

(MAE) 

054 

0.1573 

0.1380 

-0.0193 

0.1514 

-0.0060 

0.1228 

-0.0285 

0.1108 

-0.0465 

053 

0.4245 

0.3151 

-0.1093 

0.3155 

-0.1089 

0.2972 

-0.1273 

0.1399 

-0.2846 

052 

0.4528 

0.4024 

-0.0505 

0.4631 

0.0103 

0.4574 

0.0046 

0.1697 

-0.2832 

033 

0.6684 

0.4168 

-0.2515 

0.4330 

-0.2354 

0.4289 

-0.2395 

0.3431 

-0.3253 

023 

0.4358 

0.3901 

-0.0457 

0.4048 

-0.0310 

0.4048 

-0.0309 

0.3234 

-0.1123 

005 

0.5378 

0.5052 

-0.0326 

0.5409 

0.0032 

0.5413 

0.0035 

0.3928 

-0.1450 

003 

0.4909 

0.4550 

-0.0359 

0.4836 

-0.0073 

0.4837 

-0.0072 

0.3593 

-0.1316 

001 

0.6184 

0.4064 

-0.2119 

0.4644 

-0.1539 

0.4647 

-0.1537 

0.3356 

-0.2828 

Average  3D  velocity  values  and  MAE  errors  for  the  three  river  influx  test  cases 
(3DCR,  3DNR,  and  3DRR)  are  computed.  The  statistical  values  at  8  NOAA  stations  are 
recorded  in  Table  5.6.  When  compared  with  the  2D  velocity  computations,  velocities  at 
station  033  located  in  the  narrow  section  of  the  upper  bay  shows  the  greatest  improvement 
with  a  reduction  in  error  to  -0.14mA  from  -0.25mA.  This  result  concurs  with  the  previous  M2 
velocity  solution  in  which  the  3D  model  exhibits  the  most  significant  improvements  in  the 
estuarine  region  where  complex  stratification  and  mixing  due  to  tides  and  river  flux  are 
present.  The  3D  velocities  for  which  a  realistic  river  discharge  is  applied  (3DRR)  has 
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reduced  errors  in  5  of  8  locations  when  compared  to  velocitie  from  the  constant  river  flux  test 
case,  3DCR.  Time  series  of  the  velocity  magnitude  for  the  3DCR  and  3DRR  test  cases  at 
three  stations  are  shown  in  figures  5.19  through  5.21. 


Table  5.6.  Average  3D  velocities  for  modeled  and  observed  currents  and  mean  algebraic 


error  (] 

MAE)  at  8  Is 

[OAA  current  meter  sites  for  all  3D  test  cases  (3DCR,  3DNR,  and  3DRR). 

Stn 

Vavg 

(m/s) 

3DCR 

3DCR 

(MAE) 

3DRR 

3DRR 

(MAE) 

3DNR 

3DNR 

(MAE) 

054 

0.1573 

0.2216 

0.0642 

0.1374 

-0.0199 

0.0914 

-0.0660 

053 

0.4245 

0.2858 

-0.1387 

0.2777 

-0.1467 

0.2054 

-0.2191 

052 

0.4528 

0.3149 

-0.1380 

0.4722 

0.0194 

0.3593 

-0.0936 

033 

0.6684 

0.5271 

-0.1412 

0.5545 

-0.1139 

0.5275 

-0.1408 

023 

0.4358 

0.4906 

0.0549 

0.4885 

0.0527 

0.4925 

0.0567 

005 

0.5378 

0.5449 

0.0071 

0.6006 

0.0628 

0.5462 

0.0084 

003 

0.4909 

0.5283 

0.0374 

0.5410 

0.0501 

0.5297 

0.0388 

001 

0.6184 

0.4111 

-0.2072 

0.4672 

-0.1511 

0.4120 

0.2063 
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Time  series  veloctty  magnitude  at  NOS  station 005(3DCR) 


Time  KtKt  velocity  magnitude  at  NOS  station  005(3DRR) 


Figure  5.19:  Computed  (black/blue)  and  observed  (red)  current  magnitudes  at  station  005  for 
the  (a)  2DCR,  (b)  3DCR,  (c)  2DRR  and  (d)  3DRR  experiments. 
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Time  series  velocity  magnitude  at  NOS  station 0J3 


Time  series  velocity  magnitude  at  NOS  station 033(3DCR) 
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Figure  5.20:  Computed  (black/blue)  and  observed  (red)  current  magnitudes  at  station  033  for 
the  (a)  2DCR,  (b)  3DCR,  (c)  2DRR  and  (d)  3DRR  experiments. 
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Tunc  series  velocity  magnitude  at  NOS  station 054(JDCR) 


Time  series  velocity  magnitude  at  NOS  station 054 


Figure  5.21 :  Computed  (black/blue)  and  observed  (red)  current  magnitudes  at  station  054  for 
the  (a)  2DCR,  (b)  3DCR,  (c)  2DRR  and  (d)  3DRR  experiments. 
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h.  Summary 


From  the  skill  assessment  of  the  2D  and  3D  model  experiments  performed  in  this 
section,  the  findings  are  summarized  as  follows: 

1)  Water  Level  and  Current  Skill:  Given  properly  specified  boundary  conditions,  the 
ADCIRC  model  is  able  to  simulate  water  level  and  currents  reasonably  well  in  a 
tidally  driven  estuary.  The  2D  modeled  water  level  amplitudes  have  excellent 
agreement  and  indicate  little  sensitivity  to  the  specified  river  discharge  condition.  The 
3D  modeled  elevations  also  have  a  high  agreement  index  with  the  exception  of  the 
nonlinear  quarter  diurnal  constituents,  M4  and  MN4.  Both  2D  and  3D  solutions  have 
excellent  phase  propagation  characteristics. 

2)  Impact  of  River  Discharge:  The  modeled  solutions  are  quite  similar  whether  a 
constant  discharge  or  no  river  discharge  is  specified.  A  time-varying  realistic  river 
discharge  leads  to  better  model  water  level  predictions,  particularly  in  the  estuarine 
region  of  the  bay.  Velocity  computations  are  all  within  5%  error  of  one  another 
regardless  of  the  river  discharge  specification. 

3)  Influence  of  Wind  Stress:  The  application  of  wind  stress  strongly  impacts  the 
modeled  water  level  and  currents.  A  typical  strong  wind  event  of  12mA  produces  a  - 
25 cm  to  25 cm  variation  when  compared  with  a  no  wind,  tidal-driven  solution.  The 
magnitude  of  the  deviation  depends  on  the  wind  speed  and  direction. 

4)  Sensitivity  of  Bottom  Friction  Coefficient:  Despite  a  recommendation  in  the  literature 
to  increase  the  bottom  friction  coefficient  along  the  axis  of  the  bay  from  offshore 
through  the  estuarine  waters  to  the  head  of  the  river,  a  significant  reduction  in  the 
magnitudes  of  the  elevation  and  currents  led  to  poor  comparison  with  observations. 
Note  that  the  form  of  the  turbulence  closure  in  the  3D  model  was  not  modified  from 
the  constant  eddy  viscosity  model.  Changes  in  the  form  of  the  vertical  mixing  could 
significantly  alter  the  outcome.  The  interplay  between  bottom  friction  and  vertical 
mixing  was  not  explored.  The  model  sensitivity  to  the  specified  bottom  stress  for  the 
model  configurations  considered  here  is  quite  acute. 

5)  2D  vs. 3D  Solutions:  While  the  2D  model  adequately  simulates  water  level  and 
velocity  in  the  open  waters  outside  the  bay  and  in  the  river  dominated  region,  the  3D 
model  is  required  to  properly  simulate  the  complex  stratification/mixing  processes 
due  to  tidal-river  flux  interactions.  As  such  the  3D  model  computed  water  level  and 
currents  simulation  for  the  case  of  realistic  river  discharge  out  performed  the  2D 
computations.  In  the  cases  of  no  river  or  a  constant  river  discharge,  the  3D  model 
solution  did  not  provide  any  additional  accuracy.  The  simplistic  constant  eddy 
viscosity  model  used  for  vertical  mixing  could  have  reduced  the  impact  of  the  3D 
computations.  Alternate  closure  formulations  while  available  in  ADCIRC  were  not 
explored.  No  spatial  bias  in  the  error  patterns  is  discerned.  The  largest  errors 
(identified  with  three  separate  stations)  are  likely  due  to  under-resolution  or 
misrepresentation  of  small-scale  geometric  features  in  the  computational  mesh 
associated  with  those  individual  stations. 
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6.  Validation  Test  Results  -  Rattray  Island  Tidal  Currents 


Tidal  flow  around  Rattray  Island  (30°00’,S,  148°38’E)  in  the  Great  Barrier  Reef, 
Australia  (figure  6.1a)  is  an  ideal  test  case  for  the  following  reasons:  1)  an  extensive  field 
survey  has  obtained  in-situ  measurements  of  velocity  and  elevation  (Wolanski  et  ah,  1984), 
2)  the  island  lies  in  well-mixed  water  so  a  barotropic  model  is  appropriate  and  yet,  3)  the 
flow  is  strongly  three-dimensional.  Rattray  Island  is  1500m  long  and  300m  wide  and  lies  in 
water  approximately  25m  depth.  Rattray  Island  lies  perpendicular  to  tidal  currents  and  stable 
eddies  develop  in  the  wake  of  rising  and  falling  tides.  The  circulation  dynamics  at  Rattray 
Island  feature  flow  separation  and  recirculation  with  the  formation  of  eddies.  Turbid  water  in 
the  wake  of  Rattray  Island  suggest  sediment-laden  water  is  carried  upwards  to  the  surface  by 
vertical  transport  during  the  life  span  of  the  eddies  (Figure  6.1ba). 


Figure  6.1  (a)  Map  locating  Rattray  Island  in  the  Great  Barrier  Reef,  AU  (from  White  et  ah, 
2008)  and  b)  an  aerial  photo  of  Rattray  Island  taken  from  the  east  at  rising  tide  (from  Blaise 
and  White,  2006). 

Twenty-four  Aanderaa  current  meters  were  deployed  along  4  transects  in  the  wake  of 
Rattray  Island.  Velocities  were  recorded  over  the  water  column  every  10  minutes  in  the  wake 
of  Rattray  Island  (circles  in  figure  6.2)  during  a  rising  tide  in  December  of  1982.  Data  from 
current  meter  15  is  missing  and  two  of  the  current  meters,  25  and  26,  recorded  continuous  1- 
min  averaged  current  data  for  two  half-days  only.  The  remaining  current  meters  cover  a 
period  of  14  days.  Sea  level  measurements  from  three  tide  gauges  set  in  10m  of  water  around 
the  island  were  recorded  for  eight  days  starting  Nov  23,  1982  (squares  in  figure  6.2). 

From  Wolanski  et  al.  (1984),  a  3.5m  spring  tide  is  experienced  at  Rattray  Island  and 
the  rising  tide  flows  to  the  southeast  in  figure  6.1a  filling  the  bay.  A  large  time  lag,  up  to  two 
hours,  exists  between  slack  water  and  current  reversals.  The  sea  level  difference  across 
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Rattray  Island  is  phase  locked  with  tidal  currents  and  is  up  to  2.5 cm  in  amplitude.  Sea  levels 
are  higher  on  the  upstream  side  than  the  downstream  side. 

Ambient  currents  follow  the  sea  level  fluctuations  closely  though  lagging  by 
approximately  two  hours.  The  ambient  tidal  current  reverse  direction  at  an  angle  of  180°, 
making  the  ellipses  rectilinear.  The  largest  and  smallest  ranges  of  the  tidal  fluctuations  are 
3.1  m  and  1.2m  corresponding  to  the  maximum  flood  currents  at  meter  2  of  0.6mA  and 
0.4mA,  respectively.  These  currents  also  indicate  a  weak  anti-clockwise  rotation  of  the 
current,  leading  to  a  residual  circulation. 

Currents  are  deflected  around  the  island  and  show  the  generation  of  an  eddy  near  the 
lee  of  the  island  that  grows  progressively  in  size.  Typically,  after  1  hour  into  the  flood  current 
cycle,  the  current  direction  close  to  and  in  the  lee  of  the  island  (stations  24-26)  starts  to  rotate 
clockwise  until  it  reached  a  fixed  direction,  largely  against  the  ambient  current.  Progressively 
(first  at  stations  5-7  and  then  stations  12-14)  a  similar  effect  is  felt  downstream.  Thus  a 
clockwise  circulation  is  established  in  the  lee  of  the  island.  The  current  structure  in  the  eddy 
is  stable  with  time.  The  width  (parallel  to  the  long  axis  of  the  island)  of  the  eddy  along  the 
section  of  stations  1-8  is  roughly  equivalent  to  that  of  the  island  and  independent  of  the 
amplitude  of  the  ambient  current.  The  eddy  at  most  extends  to  the  section  of  stations  14-21 
yielding  a  maximum  eddy  length  of  twice  the  length  of  the  island. 

As  for  the  three-dimensional  structure  of  the  currents,  observations  near  station  24 
indicate  no  marked  velocity  gradients  over  the  water  column.  However,  a  clockwise  rotation 
of  up  to  30°  prevailed  from  surface  to  bottom. 


Figure  6.2:  A  map  of  tidal  gauge  and  current  meter  locations  deployed  in  December  1982  in 
the  wake  of  Rattray  Island  (from  Wolanski  et  al.,  1984). 
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a.  Model  Configuration 


As  part  of  the  Rattray  Island  benchmark  test  case,  five  resolutions  of  unstructured 
triangular  meshes  are  provided  and  detailed  in  Table  6.1.  The  domain  has  been  rotated  from  a 
northwest  to  southwest  orientation  to  a  north-south  orientation  so  as  to  minimize  the  x- 
component  of  the  far-field  velocity  which  will  be  used  as  a  boundary  condition.  The 
bathymetry  data  included  as  part  of  the  benchmark  has  been  interpolated  to  each  mesh.  Water 
depths  are  quite  shallow  ranging  from  9m  to  41  m  with  an  average  depth  of  25 m  (see  figure 
6.3). 


TABLE  6.1  Specifications  of  the  Rattray  Island  Benchmark  finite  element  meshes. 


Mesh  Name 

Mesh  Parameters 

Resolution 

No.  Nodes 

No.  Elements 

Min 

Max 

01 

1554 

3024 

140 

921 

02 

2531 

4973 

116 

944 

03 

3098 

6096 

88 

944 

04 

3845 

7555 

62 

721 

05 

7084 

14024 

33 

770 
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Rattray  Island  Mesh  04  With  Bathmetry 
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Figure  6.3:  Bathymetry  (in  m)  and  the  unstructured  finite  element  mesh  (04)  domain  for  the 
Rattray  Island  benchmark. 

The  only  applied  forcing  comes  from  tidal  variations  at  the  open  boundaries  (referred 
to  as  the  lower  and  upper  boundaries  (along  the  lines  y  =  -2000  m  and  y  =  12,300  m).  Tidal 
potential  forces  are  ignored  for  such  a  small  domain  and  no  surface  wind  stress  is  applied. 
Since  tidal  ellipses  are  strongly  rectilinear  in  a  direction  perpendicular  to  the  elongated  sides 
of  the  island,  closed  conditions  (IBTYPE  =  1)  are  assigned  along  the  lateral  boundaries 
(x  =  0  m  and  v  =  8000  m)  following  the  recommendations  of  the  benchmark.  For  the  upper  and 
lower  open  water  boundaries,  the  benchmark  suggests  using  the  supplied  tidal  elevation 
measured  close  to  the  island  in  combination  with  the  measured  velocities  from  a  far  field 
current  meter  station  (i.e.  23).  The  provided  elevation  and  velocity  data  start  at  day  5.11804 
and  continue  for  about  8  days.  Since  the  phase  lag  between  the  upper  and  lower  boundaries  is 
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less  than  20  min  (White  et  al.  2008),  uniform  conditions  are  applied  at  both  open  water 
boundaries. 

A  number  of  boundary  conditions  are  applied  and  evaluated  to  maximize  model 
fidelity  as  compared  to  the  observed  tidal  elevations  and  current  magnitudes.  The  results  are 
discussed  below.  The  first  implementation  of  the  open  water  boundary  forcing  involved  the 
use  of  ADCIRC’s  non-periodic  elevation  boundary  forcing,  a  fort.  19  file,  in  which  the 
provided  hourly  elevation  time  series  was  applied  at  the  lower  and  upper  boundary  nodes. 
Elevations  produced  by  this  model  reproduced  surface  elevation  magnitudes  fairly  well,  but 
produced  unrealistically  small  currents  (<  10 cm/s  maximum).  For  a  second  implementation 
of  the  model  boundary  forcing  a  tidal  boundary  flux  was  created  from  a  combination  of  both 
the  elevation  and  current  U  and  V  component  time  series  data  provided.  Model  computations 
using  this  flux  forcing  produced  surface  elevation  time  series  in  which  elevation  magnitudes 
exceeded  the  magnitude  of  the  original  tidal  elevations  (ranging  from  0  to  about  3  m)  by  a 
factor  of  3  regardless  of  the  other  model  input  parameters  such  as  vertical  eddy  viscosity 
coefficients.  Subsequent  experiments  modified  the  tidal  elevation  time  series  used  in  the  flux 
computation  by  subtracting  the  mean  from  the  entire  time  series  to  produce  a  residual 
elevation  data  set.  The  resulting  boundary  flux  forcing  data  set  is  hereafter  referred  to  as  the 
“residual”  flux.  The  modeled  surface  elevations  were  reduced  significantly  without  an 
accompanying  reduction  in  current  speed.  A  related  boundary  flux  was  computed  using  half 
of  the  values  of  the  “residual”  elevations  along  with  provided  current  data,  named  “residual- 
by-2”.  The  modeled  elevations  resulting  from  an  application  of  this  latest  boundary  flux 
forcing  produced  an  even  greater  reduction  in  the  range  of  surface  elevations,  but  also 
reduced  the  maximum  current  speed.  To  all  boundary  flux  forcing,  a  12-day  linear  or 
hyperbolic  tangent  ramping  function  was  applied,  prior  to  application  in  the  model.  As  a 
result,  the  model’s  internal  hyperbolic  ramp  function  was  not  applied.  This  latter  (“residual- 
by-2”)  and  most  successful  boundary  flux  forcing  was  assigned  through  the  fort.20  file  to  the 
northern  and  southern  boundaries  using  the  boundary  type,  IBTYPE  =  22. 

Damping  of  the  velocities  when  using  a  flux-specified  boundary  condition  have  been 
observed  in  other  ADCIRC  applications,  particularly  if  the  solution  reflects  off  the  open 
ocean  boundary  and  energy  is  unable  to  radiate  out  of  the  domain.  In  later  versions  of 
ADCIRC,  a  flux-specified,  radiative  boundary  condition  is  available  to  eliminate  this  issue 
(IBYTPE  =  52).  Sensitivity  tests  on  the  lateral  boundary  specification  indicated  that  in 
changing  from  an  impermeable  (IBTYPE  =  0)  to  a  radiative  boundary  (IBTYPE  =  30),  the 
computed  elevations  experience  a  dramatic  decrease  in  maximum  elevation  range,  from  9.1m 
to  2.1m.  Secondly,  in  comparing  the  model  response  to  flux  forcing  specified  through  an 
essential  (IBTYPE  =  2)  or  a  natural  boundary  condition  (IBTYPE  =  22),  the  essential 
boundary  condition  was  found  to  be  more  restrictive  leading  to  a  greater  potential  for  model 
instabilities.  While  highly  advective  flows  can  be  strongly  influenced  by  the  selection  of  the 
GWCE  weighting  parameter,  TAU0,  more  recent  versions  of  the  ADCIRC  model  allow  for 
spatially  varying  TAU0  values  to  be  computed  internally  in  the  model  based  on  local  velocity 
values.  However,  it  was  determined  that  the  availability  of  this  feature  would  not  lead  to 
significant  alterations  of  the  Rattray  Island  computed  elevations  and  currents. 

The  ADCIRC  model  fort.  15  input  file  parameter  specifications  for  the  3D  barotropic 
with  flux  forcing  is  configured  using  the  Makefl5  GUI.  The  default  settings  in  the  Makefl5 
GUI  were  changed  in  the  following  manner: 
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■  Activate  the  non- fatal  error  override,  NFOVER  =  1 

■  The  coordinate  system  is  Cartesian,  ICS  =  1 

■  The  model  is  run  in  3D  barotropic  mode,  IM  =  1 

■  Wetting  and  drying  is  inactive  but  the  nonlinear  finite  amplitude  terms  are  on, 
NOLIFA  =  1 

■  All  advective  terms  are  activates,  NOLICA  =  1  and  NOLICAT  =  1 

■  Specify  a  spatially  constant  Coriolis  parameter,  NCOR  =  0 

■  No  meteorological  forcing  is  applied,  NWS  =  0 

■  No  ramp  function  is  used,  NRAMP  =  0  (ramping  of  the  forcing  is  handled 
externally) 

■  A  smaller  GWCE  weighting  parameter  is  specified,  TAUO  =  0.005 

■  The  length  of  simulation  is  set  to  20  days,  RNDAY  =  20.0 

■  The  ramp  period  set  to  a  duration  of  0  days,  DRAMP  =  0.0 

■  The  minimum  depth,  HO,  is  set  to  3.0m,  HO  =  3.0 

■  The  central  projection  points  are  specified  for  the  model  domain, 
SLAM0,SFEA0  =  -88.5,  29.0,  respectively. 

■  Parameters  for  the  nonlinear  quadratic  bottom  friction  coefficient  are 
specified,  CF  =  0.05 

■  The  lateral  eddy  viscosity  coefficient  is  set  to  5,  ESL  =  5.0 

■  A  constant  Coriolis  factor  is  computed,  CORI  =  0.0000498665369 

■  Select  station  elevation  output,  NOUTE  =  1  for  just  under  8  days  (TOUTSE  = 
12.0,  TOUFE  =  20.0  model  days)  every  10  minutes  (NSPOOLE  =  300) 

■  Identify  26  elevation  station  locations,  NSTAE  =  26 

■  Select  station  velocity  output,  NOUTY  =  1  for  just  under  8  days  (TOUTSY  = 
12.0,  TOUFV  =  20.0  model  days)  every  10  minutes  (NSPOOLV  =  300) 

■  Identify  26  velocity  station  locations,  NSTAV  =  26 

■  Select  global  elevation  output,  NOUTGE  =  1  for  just  under  8  days 
(TOUTSGE  =  12.0,  TOUFGE  =  20.0  model  days)  every  10  minutes 
(NSPOOLGE  =  300) 

■  Select  global  velocity  output,  NOUTGV  =  1  for  just  under  8  days  (TOUTSGV 
=  12.0,  TOUFGV  =  20.0  model  days)  every  10  minutes  (NSPOOLGV  =  300) 

■  Quadratic  slip  is  selected,  ISLIP  =  2 

■  Surface  and  bottom  roughness  are  specified,  Z0S  =  0.01,  Z0B  =  0.005 

■  A  uniform  vertical  grid  with  21  levels  is  specified,  IGC  =  1,  NFEN  =  21 

■  A  vertical  eddy  viscosity  code  for  a  constant  is  selected,  IEVC  =  1,  and  the 
value  is  EVCON  =  0.001  (the  minimum  EV  value  is  not  used  for  IEVC  =  1) 

■  Select  3D  station  velocity  output,  I3DSV  =  1  for  just  under  8  days  (T03DSVS 
=  12.0,  T03DSVF  =  20.0  model  days)  every  10  minutes  (NSP03DSV  =  300) 

■  Identify  26  velocity  station  locations,  NSTA3DV  =  26 

■  Select  3D  global  velocity  output,  I3DGV  =  1  for  8  days  (T03DGVS  = 
1 1.9994,  T03DGVF  =  20.0  model  days)  every  hour  (NSP03DGV  =  1800) 

To  summarize,  a  fully  nonlinear  3D  model  simulation  extends  for  12  days  during  which  the 
externally  ramped  forcing  flux  is  applied  to  the  lower  and  upper  boundaries.  The  simulation 
then  continues  for  an  8  day  period  that  coincides  with  elevation  and  current  meter 
observations.  The  horizontal  mesh  is  the  highly  refined  04  mesh  whose  resolution  ranges 
from  62 m  to  721m.  The  vertical  grid  is  composed  of  21  uniform  vertical  levels.  The 
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turbulence  closure  is  a  simple  constant  vertical  eddy  viscosity  mixing  coefficient  with  a  value 
of  0.001  and  quadratic  slip  conditions  apply  at  the  sea  bed.  A  time  step  of  2  seconds  and  a 
minimum  depth  of  3  meters  is  specified. 

b.  Predictor-Corrector  Time  Stepping  Performance 

The  Rattray  Island  benchmark  provided  an  opportunity  to  test  the  predictor-corrector 
(P-C)  time  stepping  algorithm  for  a  non-wetting  and  drying  problem.  Note  applications  with 
wetting  and  drying  are  limited  by  the  propagation  speed  of  the  wetting  front.  The  wetting 
front  cannot  advance  more  than  a  single  element  per  time  step.  As  a  consequence,  the 
predictor-corrector  time  stepping  is  not  usually  an  appropriate  choice  under  wetting  and 
drying  conditions. 

The  Rattray  Island  benchmark  was  run  in  a  2D  mode  for  tests  of  the  predictor- 
corrector  time  stepping.  Recall  from  Section  1  that  the  predictor-corrector  algorithm  allows 
larger  time  steps  to  be  taken  by  the  model  for  a  give  application  and  is  activated  by  setting  a 
negative  value  for  the  time  step.  For  Rattray  Island,  the  non  P-C  application  has  a  time  step 
of  2.0  seconds.  The  model  was  run  for  successive  multiples  of  the  time  step  (i.e.  As,  8,s',  1 2,v, 
16^,  205).  The  error  analysis  is  conducted  over  1152  time  steps  so  that  in  total  there  are 
281,088  data  values  included  in  the  error  analysis.  The  percent  error  for  elevation  is 
computed  as  the  absolute  value  of  the  non  P-C  elevation  minus  the  P-C  computed  elevation 
at  244  evenly  distributed  station  locations  (figure  6.4)  throughout  the  mesh.  No  percent  error 
is  calculated  for  velocity  values  since  the  velocity  components  are  so  small.  Results  are 
recorded  in  Table  6.2.  Note  that  model  computations  using  P-C  with  time  steps  above  12s 
became  unstable.  The  elevation  differences  in  Table  6.2  between  P-C  and  non-P-C  are  less 
than  10%  for  90%  of  the  1152  time  steps  contained  in  the  analyses  at  244  stations.  So  for  the 
Rattray  Island  benchmark,  the  P-C  time  stepping  option  allows  an  increase  in  time  step  by  a 
factor  of  6  without  degradation  of  the  computed  solution. 

TABLE  6.2:  Errors  associated  with  the  predictor-corrector  time  stepping  for  the 

Rattray  Island  benchmark. 


Percent  Elevation  Differences:  (No  P-C  minus  P-C) 


%  Elevation 
Error 

V  %  Current  Error  U 

Increase 

Delt 

Min 

Max 

Min 

Max 

Min 

Max 

2x 

4s 

2.1e-07 

19.2 

0 

0.0020 

4.7e-09 

0.01 

4x 

8s 

8.3e-07 

17.7 

0 

0.0023 

1.6e-10 

0.01 

6x 

12s 

2.6e-07 

20.7 

0 

0.0026 

3.6e-10 

0.01 
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Figure  6.4:  The  regular  grid  of  244  stations  at  which  modeled  elevation  and  velocity 
solutions  from  P-C  and  non  P-C  time  stepping  are  compared. 
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c.  3D  Currents 


The  modeled  flow  pattern  on  28  November  1982  is  shown  in  figure  6.5  at  11  hrs, 
13.3 hrs,  17. 3 hrs,  and  19.5 into  the  simulation.  At  the  rising  tide  (figure  6.5a)  flow 
towards  the  north  side  of  the  island  separates  and  flows  around  the  tips  gaining  strength  to 
approximately  0.2 m/s.  A  stable  eddy  then  forms  in  the  southern  lee  of  the  island  with 
clockwise  rotation.  A  second  eddy  forms  to  the  southwest  the  island  rotating  counter¬ 
clockwise.  Asymmetry  of  the  eddies  is  due  to  the  bathymetric  variation.  The  magnitude  of 
the  currents  is  similar  to  the  free  stream  velocity  and  ranges  from  0.2-0. 1  m/s.  Just  over  two 
hours  later  (figure  6.5b)  the  eddy  has  grown  in  size  to  twice  the  width  of  the  island.  The  long 
lag  time  reported  by  Wolanski  et  al  (1984)  is  confirmed  here  and  4  hours  later  at  17.3 hrs  into 
the  simulation  the  tide  has  reversed  and  eddy  formation  is  underway  in  the  wake  of  the 
island,  now  to  north.  A  similar  two  eddy  system  of  counter-rotating  circulation  is  set-up. 
Current  magnitudes  are  somewhat  reduced  for  the  falling  tide  over  that  of  the  rising  tide.  The 
modeled  currents  do  not  seem  to  reach  the  maximums  of  0 .6m/s  and  OAm/s  as  reported  by 
Wolanski  et  al  (1984). 

Vertical  profiles  of  the  currents  at  all  current  meter  locations  are  shown  at  the  same 
times,  11  hrs,  13.3 hrs,  \1 .3hrs,  and  19.5 hrs  in  to  the  simulation  on  28  November  1982  in 
figure  6.6.  At  station  locations  in  and  around  eddy  formation  on  the  rising  tide  (figure  6.6a 
and  figure  6.6b),  stations  24,  25,  and  26  and  stations  4,  5,  and  6,  the  slight  velocity  gradients 
disappear  as  the  eddy  grows  in  size.  At  stations  in  the  far  field  away  from  the  immediate 
eddy  formation,  e.g.,  stations  11,  12,  17,  and  18,  a  boundary  layer  is  present  in  the  velocity 
profile.  Currents  on  the  non-wake  side  of  the  island  are  uniform  over  the  water  column 
(figure  6.6d). 

Currents,  both  modeled  and  observed,  located  5m  below  the  surface  are  compared  at 
the  current  meter  stations  on  28  November  1982  at  18.2 hrs  and  20.5 hrs  and  on  29  November 
1982  at  0.5 hrs,  and  2.6hrs  in  figure  6.7.  The  modeled  currents  agree  well  with  observations 
in  direction  and  magnitude  at  the  rising  tide  (figure  6.7a)  but  the  model  begins  the  tide 
reversal  earlier  than  that  measured  (figure  6.7b).  Currents  from  the  model  on  the  non-wake 
side  of  the  island  are  shifted  in  direction  from  the  measured  current  particularly  at  locations 
away  from  the  island  (figure  6.7c).  Once  again  we  observe  the  model  leading  the  phase  on 
the  tide  reversal  (figure  6.7d). 

For  completeness  another  tidal  cycle  of  horizontal  currents,  corresponding  vertical 
profiles  and  model-data  comparisons  at  current  meter  locations  are  included  in  figure  6.8, 
figure  6.9  and  figure  6.10,  respectively,  for  04  December  1982  at  3.6  hrs,  5.6  hrs,  7.6  hrs,  9.6 
hrs,  11.6 hrs,  and  13 .6hrs. 
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Figure  6.5:  Computed  currents  on  28  December  1982  at  (a)  11  hrs,  (b)  13.3 hrs,  (c)  \13hrs, 
and  (d)  19.5 hrs. 


82 


Rattray  Island  Mesh  4:  Time  Step  66,  Time  12.4583  days 
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Rattray  Island  Mesh  4:  Time  Step  80,  Time  12.5556  days 
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Figure  6.6  Vertical  profiles  of  the  modeled  currents  at  each  of  the  26  current  meter  locations 
on  28  November  1982,  (a)  1 1  hrs  and  (b)  \2>3hrs. 
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Figure  6.6  Vertical  profiles  of  the  modeled  currents  at  each  of  the  26  current  meter  locations 
on  28  November  1982,  (c)  \1.3hrs  and  (d)  19.5 hrs. 
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5000 


Rattray  Island:  Current  vectors  for  time  332.7569  days  (1982) 


Rattray  Island:  Current  vectors  for  time  332.B542  days  (1902) 


Figure  6.7:  Modeled  (black)  and  observed  (blue)  currents  at  26  current  meter  stations  on  28- 
29  November  1982  at  (a)  18.2 hrs,  (b)  20.5 hrs,  (c)  0.5 hrs,  and  (d)  2 ,6hrs. 
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Figure  6.8:  Computed  currents  on  04  December  1982  at  (a)  3.6  hrs,  (b)  5.6  hrs,  (c)  7.6  hrs, 
(d)  9.6  hrs,  (e)  11 .6hrs,  and  (f)  13 ,6hrs. 
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Rattray  Island  Mesh  4:  Time  Step  843,  Time  17.8542  days 
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Rattray  Island  Mesh  4:  Time  Step  855,  Time  17.9375  days 
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Figure  6.9:  Vertical  profiles  of  the  modeled  currents  at  each  of  the  26  current  meter  locations 
on  04  December  1982,  (a)  3.6  hrs  and  (b)  5.6  hrs. 
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Rattray  Island  Mesh  4:  Time  Step  867,  Time  18.0208  days 
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Rattray  Island  Mesh  4:  Time  Step  879,  Time  18.1042  days 
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Figure  6.9:  Vertical  profiles  of  the  modeled  currents  at  each  of  the  26  current  meter  locations 
on  04  December  1982,  (c)  7.6  hrs  and  (d)  9.6  hrs. 


Rattray  Island  Mesh  4:  Time  Step  891,  Time  18.1875  days 
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Figure  6.9:  Vertical  profiles  of  the  modeled  currents  at  each  of  the  26  current  meter  locations 
on  04  December  1982,  (e)  11 ,6hrs,  and  (f)  13 ,6hrs. 
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Rattray  Island:  Current  vectors  for  time  338.1528  days  (1382) 


Rattray  Island:  Current  vectors  for  time  338.3194  days  (1982) 


Rattray  Island:  Current  vectors  for  time  330.4861  days  (1902) 


Rattray  Island:  Current  vectors  for  time  338.2361  days(igsz) 


Rattray  Island:  Current  vectors  for  time  336.5694  days  (1962) 


Figure  6.10:  Modeled  (black)  and  observed  (blue)  currents  at  26  current  meter  stations  on  04 
December  1982  at  (a)  3.6  hrs,  (b)  5.6  hrs,  (c)  7.6  hrs,  (d)  9.6  hrs,  (e)  11 ,6hrs,  and  (f)  13 ,6hrs. 


90 


d.  Summary 

The  application  of  ADCIRC  to  the  Rattray  Island  benchmark  demonstrates  that  the 
model  has  skill  not  only  in  the  replication  of  asymmetric  eddy  formation  in  the  wake  of  an 
island  on  the  rising  and  falling  tide  but  also  the  expected  vertical  structure  of  the  currents. 
The  computed  magnitudes  of  the  tide  and  tidal  current  were  not  able  to  match  the 
observations.  This  is  largely  thought  to  be  an  artifact  of  the  boundary  condition  specification 
of  the  benchmark  itself  and  especially  in  the  context  boundary  types  available  within  the 
ADCIRC  model. 

To  summarize  the  Rattray  Island  experience:  1)  the  selection  of  boundary  forcing  data 
and  boundary  type  is  particularly  important  for  small  domains,  and  2)  boundary  effects  can 
be  minimized  by  pushing  open  boundaries  away  from  primary  circulation  features  in  such  a 
way  that  energy  can  radiate  away  from  the  localized  region  of  interest. 

While  some  experimentation  was  conducted  with  regard  to  the  form  of  the  3D 
turbulence  closure  and  the  associated  parameter  values,  a  rigorous  study  was  deemed  beyond 
the  scope  of  the  present  validation  testing.  It  is  known  however  that  the  model  results  are 
sensitive  to  the  bottom  friction  parameterization  and  the  turbulence  closure  (vertical  eddy 
viscosity).  Recently  a  significant  level  of  effort  has  been  applied  to  improve  the  3D 
computations  in  the  model.  Recent  versions  of  the  code,  i.e.  v47  and  higher,  have  made  bug 
corrections  and  systematic  improvements  to  the  Mellow-Yamada  2.5  turbulence  closure 
implementation  within  ADCIRC  with  the  specific  goal  to  make  it  more  robust  in  shallow 
water.  For  example,  the  turbulent  length  scale  initialization  has  been  modified  and  new  limits 
defined  for  the  minimum  and  maximum  value  of  the  turbulent  length  scale.  Additional 
modifications  have  targeted  the  bottom  drag  coefficient  calculation. 

Lastly,  more  recent  versions  of  the  ADCIRC  model,  v47  and  higher,  have  available  a 
radiative  flux  specified  boundary  condition.  This  boundary  type  would  likely  lead  to  better 
model  performance  for  small  domain,  shallow  water  problems  like  Rattray  Island. 

7.  Validation  Test  Results  -  Forecasting  for  MREA07 

Sea  trials  conducted  in  the  shallow  Gulf  of  LaSpezia  in  collaboration  with  the  NATO 
Undersea  Research  Center  (NURC)  as  part  of  the  Marine  Rapid  Environmental  Assessment 
2007  (MREA07)  provided  an  opportunity  to  test  and  evaluate  all  aspects  of  the  software 
under  transition  to  NAVOCEANO  as  well  as  the  performance  of  the  model  predictions 
themselves.  The  developed  enabling  software,  MeshGUI,  Makefl5,  and  Makef22  were  all 
used  to  set  up  a  2D  ADCIRC  model  simulation  for  the  region.  Scripting  software  was  then 
created  to  run  the  ADCIRC  model  in  a  simulated  “real-time”  prediction  scenario,  producing 
marine  environmental  forecasts  for  the  two-week  period  spanning  observational  cruises  in 
the  summer  of  2007.  Below  we  describe  implementation  of  the  enabling  software  for  the 
Gulf  of  La  Spezia,  define  parameters  of  the  “real-time”  ADCIRC  forecast  simulations  and 
compare  the  model  to  observations  taken  as  a  part  of  the  coastal  experiment,  POET 
(Predictive  Oceanographic  Experimental  Trial). 

The  MREA07  focused  on  the  Gulf  of  La  Spezia  located  on  the  eastern  side  of  the 
Ligurian  coast  (northwest  Italy)  (Figure  7.1).  Delimited  by  the  Tino  and  Palmaria  islands  on 
the  west  and  by  Punta  Bianca  peninsula  on  the  east,  the  Gulf  is  surrounded  by  mountains. 
Northwest-southeast  oriented,  it  is  5  km  wide  while  its  length  is  about  10  km.  A  dam  (length 
of  2.2  km)  separates  the  Gulf  in  two  areas:  inside  the  dam  there  is  the  harbor  having  a  mean 
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depth  of  about  10-11  m;  outside  the  dam,  representing  transition  to  the  open  sea,  is  deeper 
and  has  an  irregular  bathymetry.  Depth  of  the  Gulf  progressively  increases  in  a  westward 
direction:  the  maximum  depth  (about  25  m)  is  in  proximity  of  Palmaria  and  Tino  islands, 
where  the  bathymetry  becomes  very  steep.  Two  passages,  at  the  dam  ends,  permit  exchange 
between  the  inside  and  outside  part  of  the  Gulf.  The  western  opening  is  wider  (360  m)  and 
deeper  (15-16  m  )  than  the  eastern  one  (180  m  width  and  11-12  m  depth).  While  the  main 
connection  with  the  open  sea  is  through  the  Gulf  open  boundary  and  the  passage  between 
Palmaria  and  Tino  islands,  some  exchange  is  also  possible  through  the  Portovenere  channel 
(between  land  and  Palmaria  island),  although  the  channel  is  very  narrow  (150  m  wide)  and 
very  shallow  (the  sill  depth  is  3  m). 


Figure  7. 1 .  A  map  of  the  Gulf  of  La  Spezia. 

a.  Model  Configuration 

To  model  the  Gulf  of  La  Spezia  using  ADCIRC,  a  new  unstructured  mesh  is  created 
using  the  MeshGUI  tool.  Required  inputs  to  the  MeshGUI  tool  include  an  ordered  list  of 
coastal  outline  segments  and  bathymetry.  The  coastal  outline  was  extracted  from  the  World 
Vector  Shoreline  database  (http://rimmer.ngdc.noaa.gov/)  and  provided  by  NURC  through 
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the  MREA07  collaboration.  Bathymetric  data  is  of  nominal  3  km  resolution  and  is  a  sub¬ 
sampled  composite  from  Italian  and  French  hydrographic  databases,  also  provided  by  NURC 
through  the  MREA07  collaboration  (Figure  7.2).  A  strategy  that  places  the  open  ocean 
boundary  away  from  Gulf  of  La  Spezia  coastal  waters  is  used  such  that  the  model  domain 
includes  all  of  the  offshore  and  nearshore  waters  and  embayments  circumscribed  by  the 
seaward  boundary  and  the  coastline  shown  in  figure  7.3.  Upstream  portions  of  the  Magra 
River  were  truncated  from  the  mesh  to  omit  extremely  shallow  depths  (less  than  1  m)  from 
the  modeled  region.  The  westernmost  and  southern  boundaries  are  configured  as  a  straight 
line  to  facilitate  the  extraction  of  boundary  forcing  data  from  the  nested  Navy  Coastal  Ocean 
Model  (NCOM).  Other  grid  generation  parameters  specified  within  MeshGUI  are  the 
maximum  offshore  resolution  of  200  km  and  the  number  of  refinements,  set  to  4.  The 
resulting  mesh  (figure  7.3)  has  24,740  nodes  and  47,369  triangular  elements  and  spatial 
resolution  ranging  from  16  m  to  181  m.  This  mesh  was  cast  in  the  form  of  a  fort.  14  grid  file. 
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MREA07  AS  14  Region:  Bathymetry,  Coastline  (blue  points),  Boundary  (black  circles)  Depth  (m) 
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Figure  7.2.  Bathymetry  contours  for  the  Gulf  of  La  Spezia. 
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Figure  7.3.  The  Gulf  of  La  Spezia  unstructured  mesh  created  by  MeshGUI.  Resolution  ranges 
from  16  m  to  181  m. 


Determination  of  the  appropriate  forcing  for  the  model  is  the  next  step.  Entering  into 
the  model  domain  are  three  sources  of  discharge  (figure  7.4),  daily  outflow  from  the  ENEL 
power  plant  equal  to  30  m3/s,  the  ENEL  cooling  pump  daily  discharge,  also  30  m3/s,  and  the 
daily  climatological  average  discharge  for  the  Magra  River  shown  in  figure  7.5  for  June  and 
July.  Discharge  into  the  model  domain  (specified  as  IBYTPE  22)  enters  through  a  fort.20  file 
which  is  created  outside  of  the  Makefl5  GUI.  Forcing  for  the  open  ocean  boundaries  as 
previously  mentioned  is  taken  from  the  innermost  3  km  nest  of  a  3-level  nested  NCOM 
model  of  the  Ligurian  Sea,  that  includes  La  Spezia  Bay.  The  nested  domains  of  the  NCOM 
MREA07  implementation  are  shown  in  figure  7.6.  The  NCOM  model  included  surface 
forcing  from  heat  flux,  wind  and  pressure  derived  from  the  Coupled  Ocean  Atmosphere 
Mesoscale  Prediction  System  (COAMPS),  Europe_2,  at  27  km  resolution.  Tides  for  the 
NCOM  model  are  applied  at  the  open  water  boundaries  of  the  second  nest  (9  km)  using 
values  from  the  Oregon  State  global  tidal  database  (Egbert  and  Erofeeva,  2003).  Tidal 
variability  for  the  innermost  NCOM  model  nest  is  forced  only  through  the  NCOM  solution  of 
the  second  nest.  Daily  surface  elevation  and  U  (eastward)  and  Y  (northward)  current 
components  available  at  1-hour  intervals  are  extracted  from  the  NCOM  model  from  the  24-hr 
hindcast  (ADCIRC  analysis)  and  the  48-hr  forecast  period.  NCOM  current  components 
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within  the  uppermost  40  m  were  depth-averaged  to  obtain  a  velocity  that  was  consistent  with 
the  2D  version  of  ADCIRC. 

A  normal  specified  flux  (IBTYPE  =  2)  is  applied  along  the  southern  boundary  of  the 
ADCIRC  model  domain.  Along  this  boundary,  NCOM  elevation  and  velocities  previously 
extracted  at  the  NCOM  grid  points  (shown  as  red  squares  in  figure  7.4)  are  interpolated  onto 
the  nearest  neighbor  ADCIRC  mesh  nodes  (shown  as  blue  crosses  in  figure  7.4).  Flux  values 
are  computed  using  standard  relationships  between  elevation,  velocity,  and  the  normal  flux 
through  a  segment  to  obtain  the  normal  flux  per  unit  width  in  units  of  m2/s.  Along  the 
western  boundary,  radiation  boundary  conditions  are  selected  (IBTYPE  =  30)  to  allow  non- 
modeled  energy  to  exit  the  domain  (no  actual  forcing  values  are  applied).  Figure  7.7  clearly 
delineates  the  boundary  types  specified  over  the  ADCIRC  model  domain. 


MREA07:  ASI4  Region,  GP's  NCOM3  Nest  (Squares)  Near  ADCIRC  Open  Boundary  (Points) 
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Figure  7.4.  Freshwater  discharge  locations  of  the  ENEL  Outflow,  the  ENEL  Cooling  Plant 
Pump,  and  the  Magra  River.  The  ENEL  daily  discharges  are  each  specified  as  30  m3/s. 
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Magra  River  Daily  Average  Discharge,  1951-1965 
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Figure  7.5.  Climatological  (1951-1965)  daily  average  discharge  values  for  the  Magra  River 
during  the  months  of  June  and  July. 
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MREA07:  ASI4  Region,  GP:s  NCOM  Nests,  Discharges 
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Figure  7.6.  NCOM  model  nests  of  km  (red),  km  (green)  and  km  (blue)  in  the  Ligurian  Sea. 
Elevation  and  currents  were  applied  from  the  NCOM3  nest  to  the  ADCIRC  southern 
boundary. 
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Figure  7.7.  Boundary  type  specifications  for  the  ADCIRC  model  domain.  Radiation 
conditions  (IBTYPE  =  30)  are  depicted  in  red  while  the  normal  flux  boundary  (IBTYPE  = 
30)  is  depicted  in  blue. 

Surface  meteorological  forcing  is  applied  using  the  highest  resolution  products 
available  at  the  time  of  the  MREA07,  the  ALADIN  products  made  available  to  MREA07 
collaborators  by  EPSHOM,  the  French  Navy’s  Hydrographic  and  Oceanographic 
Department.  The  ALADIN  products  for  mean  sea  level  pressure,  wind  magnitude,  and  wind 
direction  have  0.1  degree  resolution  with  the  data  arriving  in  two  parts:  the  first  part 
contained  hourly  forecasts  from  1  to  12  hr,  and  the  second  contained  three-hourly  forecasts 
from  15  to  48  hr  (figure  7.8).  Note  that  there  was  no  analysis  at  0  hr.  The  second  part  of  the 
data  set  was  interpolated  to  hourly  intervals,  and  a  “pseudo  analysis”  was  created  by 
interpolating  between  the  current  day’s  hour  1  forecast  and  the  previous  day’s  hour  23 
forecast.  In  this  manner,  a  data  set  at  each  day’s  0  hr  was  created  to  provide  continuous 
temporal  coverage.  This  continuous  time  series  of  ALADIN  wind  products  were  processed 
by  the  Makef22  software  to  create  the  fort.22  meteorological  forcing  file  for  the  ADCIRC 
model  simulation.  The  Makef22  utility  handed  the  data  read,  interpolation  of  the  data  from  a 
regular  rectangular  grid  onto  the  ADCIRC  finite  element  mesh,  and  export  of  the  interpolated 
data  to  a  fort.22  formatted  file.  The  ramping  of  the  wind  forcing  for  the  initial  15-day  ramp- 
up  period  is  also  handled  by  the  Makef22  software. 
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MREA07:  ASI4  Region,  ALADIN  Winds  lor  200705030600_01  5 


Figure  7.8.  ALADIN  wind  vectors  over  La  Spezia  Bay.  Resolution  is  0.1  degree  with  an 
hourly  frequency  over  the  first  12  hours  and  3 -hourly  frequency  out  to  48  hours. 


The  final  step  in  the  preparation  of  an  ADCIRC  model  simulation  is  the  creation  of 
the  parameter  and  periodic  forcing  file,  fort.  15.  The  Makefl5  GUI  was  applied  for  this 
purpose.  The  default  settings  in  the  Makefl5  GUI  were  changed  in  the  following  manner: 

■  Activate  the  non- fatal  error  override,  NFOVER  =  1 

■  Wetting  and  drying  and  the  finite  amplitude  term  contributions  are  disabled, 
NOLIFA  =  0 

■  Specify  a  spatially  constant  Coriolis  parameter,  NCOR  =  0 

■  Activate  tidal  potential  forcing,  NTIP  =  1 

■  Time  step  is  reduced  to  1  sec,  DTDP  =  1.0 

■  The  reference  time  is  set  to  the  time  after  ramping,  REFTIM  =15.0 

■  The  length  of  simulation  is  set  to  19  days,  RNDAY  =  17.0 

■  The  ramp  period  is  set  to  a  duration  of  15  days,  DRAMP  =  15.0 

■  The  central  projection  points  are  computed  for  the  model  domain, 
SLAM0,SFEA0 

■  A  constant  Coriolis  factor  is  computed,  CORI  =  0.0000252202771 

■  Assign  1  tidal  potential  constituent,  NTIF  =  1 

■  Tidal  potential  constituent,  TIPOTAG  =  M2 

■  Input  the  date  at  the  start  of  the  simulation  to  compute  the  nodal  factors  (June 
1,  2007) 

■  Select  global  elevation  output,  NOUTGE  =  -1  for  2  days  (TOUTSGE  =  15.0, 
TOUFGE  =  17.0  model  days)  every  hour  (NSPOOLE  =  3600) 

■  Select  global  velocity  output,  NOUTGV  =  -1  for  2  days  (TOUTSGV  =  15.0, 
TOUFGV  =  17.0  model  days)  every  hour  (NSPOOLV  =  3600) 
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Output  a  hotstart  at  the  24th  hour  of  the  forecast  period,  NHSTAR  =  1, 
NHSINC  =  1378800  (model  day  =  15.958333) 


All  of  the  needed  input  files  for  an  ADCIRC  model  simulation  are  described,  fort.  14, 
fort.  15,  fort.20,  and  fort.22.  To  summarize,  the  model  configuration  will  produce  water  levels 
and  2D  currents  subject  to  surface  wind  and  pressure  forcing,  discharge  from  two  locations  at 
the  ENEL  power  plant  as  well  as  the  Magra  River,  and  offshore  oceanic  forcing  represented 
by  the  NCOM  model  solution.  Tidal  potential  forcing  within  the  ADCIRC  domain  is  also 
included  and  Coriolis  effects  remain  constant  over  the  domain.  The  only  nonlinearity 
activated  in  the  model  is  the  bottom  stress  term.  The  predictive  forecast  system  configured 
using  these  input  files  is  outlined  in  the  next  section. 

b.  Predictive  System 

The  objective  of  the  MREA07  model  forecast  system  was  to  generate  48  hour 
forecasts  of  elevation  and  2D  currents  for  the  Gulf  of  La  Spezia  over  the  period  1 8  June  2007 
to  1  July  2007.  To  do  this  an  initial  model  run  was  conducted  to  spin-up  the  tidal  solution 
within  the  domain.  As  mentioned  in  the  previous  section,  this  run  started  on  1  June  2007 
applying  a  hyperbolic  ramp  function  that  lasted  for  15  days.  After  ramping  a  2-day  forecast 
begins.  At  the  23rd  hour  of  the  forecast  period  a  model  solution  is  written  as  a  hotstart  file 
(fort.67)  (model  day  15.9583333).  This  field  is  used  to  initialize  the  next  day  48-hour 
forecast.  Each  subsequent  forecast  run  includes  one  hour  prior  to  forecast  to  allow  the 
interpolation  and  application  of  the  analysis  ALADIN  wind  and  pressure  fields  valid  1  hour 
into  the  forecast.  This  out-of-sync  timing  of  the  meteorological  forcing  requires  that  the 
model  fields  be  output  prior  to  the  start  of  the  forecast  during  a  time  which  had  valid 
meteorological  forcing,  i.e.  1  hour  before  the  start  of  a  forecast  run.  Linear  interpolation 
between  the  23rd  hour  of  the  previous  day  forecast  and  1  hour  into  the  current  day  forecast 
result  in  a  seamless  application  of  the  meteorological  forcing.  These  2-day  +  1  hour  model 
simulations  continue  through  the  1th  of  July  2007,  the  time  at  which  NCOM  model 
availability  ceased. 

For  successive  model  forecast  runs,  the  following  parameters  in  the  fort.  15  files  were 
changed  via  a  script: 

■  The  length  of  simulation  is  incremented  by  1,  RNDAY  =  RNDAY+1 

■  Times  for  global  elevation  output  are  incremented  by  1,  TOUTSGE  = 
TOUTSGE+ 1 ,  TOUTFGE  =  TOUTFGE+ 1 

■  Times  for  global  velocity  output  are  incremented  by  1,  TOUTSGV  = 
TOUTSGV+1,  TOUTFGV  =  TOUTFGV+1 

■  Timestep  for  hotstart  file  generation  is  incremented  by  1  day  or  NHSINC  = 
NHSINC  +  86400  (no.  of  timesteps  per  day) 

The  49-hour  forecast  ADCIRC  model  simulation  was  executed  in  parallel  on  5  AMD 
dual-core  Opterons,  model  252,  having  a  single-core  speed  of  2.6  GHz.  The  wallclock  time 
for  a  single  forecast  run  was  1  hr  15  min  on  average.  This  amounts  to  0.00255  CPU  sec  per 
time  step  per  processor  for  this  24,740  node  grid  model. 
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c.  2D  Currents 


As  part  of  the  Predictive  Oceanographic  Environmental  Trial  (POET)  a  vessel 
mounted  RDI  300  kHz  ADCP  was  used  to  collect  data  along  the  water  column  in  the  external 
part  of  the  Gulf,  along  five  tracks  (figure  7.9).  The  ADCP  cell  size  was  1  m  and  the  sampling 
rate  was  15-sec.  Ten  minute  positions  along  the  ADCP  tracks  are  shown  as  red  circles  in 
figure  7.9.  No  detiding  procedure  has  been  applied  due  to  the  smallness  of  tidal  signal.  The 
cruise  began  on  at  13:23  GMT  in  the  north  western  corner.  A  zigzag  track  was  followed  to 
maximize  spatial  coverage  until  ending  nearly  3  14  hours  later  at  16:48  GMT  in  the  south  east 
comer  waters  bounded  by  Tino  Island  to  the  west  and  the  Punta  Bianca  peninsula  to  the  east. 


MREA07:  ASI4  Region,  ADCIRC  Open  Boundary,  POET  ADCPs  and  Discharges 


Figure  7.9.  Map  of  the  ADCP  tracks  from  the  POET  experiment  on  19  Jun  2007.  Red  circles 
indicate  10-minute  intervals. 


Data  drop  outs  ranged  from  30  sec  to  4  min  15  sec  with  139  data  gaps  greater  than  15 
sec.  To  perform  meaningful  model-data  comparisons  using  such  data  it  was  necessary  to  first 
perform  a  time  interpolation  on  the  data.  A  14  hour  moving  average  was  applied  to  smooth 
the  data  to  better  match  the  hourly  frequency  of  the  model  output.  Lastly,  the  current 
observations  were  depth-averaged  so  that  direct  comparisons  with  the  2D  ADCIRC  modeled 
currents  could  be  made.  This  processing  of  the  U  and  V  components  of  the  current 
observations  is  shown  in  figure  7.10  where  the  interpolated  data  is  depicted  by  red  lines  and 
the  14  hour  averaged  data  is  displayed  in  blue. 

The  magnitude  and  direction  as  well  as  U  and  V  components  of  the  ADCIRC 
currents  (red)  are  compared  to  the  observed  currents  (blue)  as  shown  in  figures  7.11.  The 
computed  differences  are  plotted  on  the  same  graphs  in  black  with  the  magnitudes  indicated 
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on  the  right-hand  scale.  The  trend  exhibited  by  the  observations  in  both  current  magnitude 
and  by  component  is  closely  followed  by  the  modeled  values.  Some  variability  of  the 
observations  is  not  captured  by  the  1-hour  frequency  of  the  model  data.  Statistical 
representation  of  the  model-data  error  is  provided  by  figure  7.12.  The  model-data  correlation 
with  respect  to  current  magnitude  is  over  0.97  while  the  U  and  V  component  correlation 
values  are  0.94  and  0.92,  respectively.  The  largest  error  clearly  enters  through  the  modeled 
current  direction  which  is  inversely  related  and  largely  un  correlated  to  observations.  From 
the  time  series  of  direction  in  figure  7.13,  most  of  the  directional  error  occurs  during  the  first 
two,  northernmost  legs  of  the  ADCP  transects. 

Comparisons  of  the  ALADIN  wind  speed,  direction  and  air  pressure  to  the  ENEA 
meteorological  station  data  (figure  7.14)  clearly  indicate  that  the  wind  directions  produced  by 
the  ALADIN  model  are  not  realistic.  The  ENEA  meterological  station  is  located  along  the 
eastern  shoreline  across  from  the  dam.  The  wind  field  error  would  be  cause  for  significant 
error  in  the  modeled  currents.  Note  also  that  the  two  northernmost  ADCP  segments  would  be 
using  the  ALADIN  winds  that  are  coming  from  the  land  while  remaining  three  southernmost 
ADCP  segments  would  be  using  ALADIN  wind  values  that  were  computed  over  the  water. 
This  is  a  consequence  of  the  rather  coarse  resolution  of  the  ALADIN  winds  for  the  shallow 
coastal  region  under  study. 
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Current  Magnitude  (rrVs)  Current  Magnitude  (rrVs) 


MREAQ7:  Depth -Averaged  POET  ADCP  Time  Series  (U  Component) 
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Figure  7.10.  Time-interpolated  (a)  U  and  (b)  V  components  of  the  current  observations  (red) 
and  the  14  hour  moving  averaged  data  (blue). 
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MREA07:  ADCIRC  and  POET  ADCPTime  Series,  U  Component 
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Figure  7.11.  Modeled  (red)  and  observed  (a)  U  and  (b)  V  components  of  the  current;  the 
difference  is  shown  in  black  (right-hand  scale). 
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Figure  7.12.  Model  vs.  observed  currents  (a)  magnitude,  (b)  U  component,  (c)  direction,  and 
(d)  V  component  of  the  current;  correlations  are  given  in  the  upper  left  of  each  panel. 
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MREA07:  ADCIRC  and  POET  ADCP  Time  Series,  Current  Speed 


Figure  7.13.  Modeled  (red)  and  observed  (a)  speed  and  (b)  direction  of  the  current;  the 
difference  is  shown  in  black  (right-hand  scale). 
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Figure  7.14.  Model  vs.  observed  winds  (a)  speed,  (b)  direction,  and  (c)  mean  sea  level 
pressure;  correlations  are  given  in  the  upper  left  of  each  panel. 


d.  Summary 

The  MREA07  exercise  provided  an  opportunity  to  demonstrate  application  of  the 
enabling  software,  MeshGUI,  Makefl5,  and  Makef22,  to  set-up  the  ADCIRC  model  for  a 
new  geographical  region.  The  finite  element  mesh  and  the  associated  grid  file,  fort.  14,  were 
created,  as  well  as  the  model  parameter  and  periodic  forcing  file,  fort.  15  and  the 
meteorological  forcing  file,  fort.22.  The  model  once  configured  was  run  in  a  forecast  mode 
for  two  weeks  in  June,  2007  with  forcing  applied  from  real-time  meteorological  model 
(ALADIN)  and  regional  oceanic  products  (NCOM).  Model-data  comparisons  for  2D  currents 
collected  from  an  ADCP  survey  lend  confidence  to  the  modeled  currents  and  highlight  their 
sensitivity  to  the  surface  wind  forcing. 
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